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(ABSTRACT) 

The objective of this study was to determine the effect of nonlinear material behavior on the 
response and failure of unnotched and notched angle-ply laminates under uniaxial 
compressive loading. The endochronic theory was chosen as the constitutive theory to model 
the AS4/3502 graphite-epoxy material system. 

Three-dimensional finite element analysis incorporating the endochronic theory was used to 
determine the stresses and strains in the laminates. An incremental/iterative initial strain 
algorithm was used in the finite element program. To increase computational efficiency, a 
180° rotational symmetry relationship was utilized and the finite element program was 
vectorized to run on a super computer. 

Laminate response was compared to experiment revealing excellent agreement for both the 
unnotched and notched angle-ply laminates. Predicted stresses in the region of the hole were 
examined and are presented, comparing linear elastic analysis to the inelastic endochronic 
theory analysis. 

A failure analysis of the unnotched and notched laminates was performed using the quadratic 
tensor polynomial. Predicted fracture loads compared well with experiment for the unnotched 
laminates, but were very conservative in comparison with experiments for the notched lami- 


nates. 
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1.0 Introduction and Literature Review 


1.1 Introduction 


Composite materials present increased flexibility to structural designers because laminates 
from these materials can be tailored to achieve desired stiffness, strength, and dimensional 
stability properties, A common class of laminates is the angle-ply laminate which consists 
of ± <p plies. These laminates often exhibit nonlinear material response under in-plane 
loading. The laminate with a hole is also common in many design applications. When loaded, 
such notched laminates have an inplane stress concentration at the hole boundary as well as 
significant interlaminar stresses at all free edges. Few researchers have studied the complex 
three-dimensional stress distribution in compression-loaded angle-ply laminates with a hole. 

The angle-ply laminate is commonly used in advanced aircraft and aerospace designs be- 
cause its stiffness can be varied over a wide range, dependent on the layup angle 9. Of 
course, in such design applications, cutouts are often a necessity. Therefore, it is important 
to know how cutouts affect the stiffness and strength of these laminates. Because composite 
material strength is lower under compression than in tension, most advanced designs are 
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compression loading critical. With these considerations, it is easy to see the importance of 
studying the response of notched angle-ply laminates under uniaxial compression loading. 

The objective of the present study is to predict and investigate the three-dimensional stress 
distribution in compression loaded angle-ply laminates with a central hole. A three- 
dimensional endochronic constitutive theory for transversely isotropic composite materials is 
used to model the nonlinear material response. The endochronic theory and 
incremental/iterative constant stress and constant strain algorithms are incorporated in a fi- 
nite element program to determine laminate stresses. Inversion symmetry is utilized to sig- 
nificantly reduce the computer storage requirements and execution time. Failure initiation 
and failure modes are predicted using the tensor polynomial and maximum stress failure cri- 
teria. The analytical results from this study are compared with existing experimental results. 


1.2 Literature Review 


This study involves a broad range pf subject areas, among them the endochronic constitutive 
theory, three-dimensional finite element analysis incorporating the endochronic theory, ma- 
terial nonlinear laminate response, and notched laminate strength prediction. A thorough 
review of the above areas would be exhausting; therefore, only the most significant subject 
areas are covered. The following review encompasses material nonlinearity in composites, 
the endochronic constitutive theory, and finite element analysis with material nonlinearity. 
For a review on the stress analysis and strength predictions of laminated plates with holes, 
see Burns, Herakovich, and Williams \ 
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1.2.1 Material Nonlinearity in Composites 


The literature is rich with papers describing nonlinear constitutive models for isotropic mate- 
rials. However, for composite materials, the constitutive theory must include anisotropic re- 
sponse. With this additional requirement of anisotropy, the choices of models are more 
limited. Because composites can be analyzed both on a microscopic scale (as a heteroge- 
neous medium) and on a macroscopic scale (as a homogeneous anisotropic medium), com- 
posite material constitutive theories are divided into micro-theories and macro-theories. 

1. 2.1.1 Microscopic Theories 

Microscopic constitutive theories use the elastic and inelastic mechanical properties of the 
fiber and matrix to formulate the anisotropic constitutive relations on a macroscopic level. 
Dvorak et al 2 and Aboudi 3 have made significant contributions to the development of micro- 
scopic constitutive theories. Many others have contributed as well. Aboudi 3 and Pindera and 
Herakovich 19 both give thorough reviews of the literature on microscopic constitutive theories 
for composites. 


1.2.1. 2 Macroscopic Theories 

In contrast to the microscopic approach, macroscopic constitutive theories treat the lamina 
as a homogeneous anisotropic medium and do not consider the individual constituent prop- 
erties. Many researchers have contributed to the development of macroscopic constitutive 
theories for composites. Only a limited number are included in this review. For a thorough 
review of the literature on macroscopic composite material constitutive theories, see Pindera 
and Herakovich 19 . 
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One of the earlier works in the area of nonlinear constitutive theories for fibrous composites 
was that of Hahn and Tsai 4 . They used complementary elastic energy density to develop a 
material nonlinear stress-strain relationship for a two-dimensional lamina. The stress-strain 
relationship is linear in uniaxial loading in the longitudinal and transverse directions, but 
nonlinear in shear. They compared their theory with off-axis coupon tests and obtained sat- 
isfactory results for boron/epoxy Narmco 5505 and graphite/epoxy 4612 Morganite II. 

Sandhu 5 developed an incremental constitutive theory using the tangential properties of the 
lamina. The strain-dependent tangential properties were used to express the strain increment 
in terms of the stress increments. Off-axis unidirectional and multi-directional laminate tests 
were performed to compare the theory to experiment. 

Hashin, Rosen, and Bagchi 8 proposed a theory which modeled the nonlinear behavior of a 
lamina in shear and transverse directions. The longitudinal (fiber direction) stress is modeled 
as linear elastic and does not affect the shear and transverse strains. The theory includes 
interaction between the transverse and shear stresses in the nonlinear range. 

An interesting approach was taken by Dvorak et al 7 to develop a continuum theory for the 
elastic-plastic response of fibrous composites. They used a microscopic scale finite element 
analysis on materials such as boron-aluminum to develop the theory. The theory agreed well 
with numerical analysis for a variety of loading programs. 


1.2.2 Endochronic Theory 


The endochronic theory has most often been applied to isotropic materials. Although the 
endochronic equations were not originally derived to model anisotropic response, the theory 
is of a general nature, allowing it to model many material types. Recently, the theory has 
been extended to anisotropic materials. 
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1. 2.2.1 Isotropic Materials 


Although there is some disagreement as to when and with whom the theoretical foundations 
for the endochronic theory were developed, Valanis 11 12 is credited with one of the earliest 
published works of the endochronic theory. He used the concepts of irreversible 
thermodynamics along with experimentally observed material behavior to develop a math- 


O fYl otirQ I fhorirt/ tkAKWi a t nr>/vnnl » TUn j l ' _ xl n , , 

mwi; wi moiM.v-viguupiaouLiijf, me iciin ciiuubfiiuiMu theory encompasses ail 


those theories in which "the state of stress at the present time is a function of the history of 
strain with respect to a time scale, which is not the absolute time scale measured by a clock, 
but a time scale which in itself is a property of the material at hand" 


In a series of papers, Valanis 13 - 18 and Valanis and Lee 17 extended the endochronic theory 
to model increasingly complex phenomena. They compared various theoretical studies with 
experimental results and found good agreement in most instances. Through these papers, 
they brought to light the numerous applications which could be accurately modeled with the 
endochronic theory. 


1.2. 2. 2 Anisotropic Materials 

Only two previous distinct research efforts have directed attention towards the development 
and application of the endochronic theory to fibrous composites. They are the works of 
Pindera and Herakovich 18 20 and Zinov'ev and Sarbaev * Pindera and Herakovich 13 20 ex- 
tended Valanis' endochronic theory to anisotropic materials and then used this theory to 
model the nonlinear response of unidirectional Celion 6000 graphite-polyimide. They derived 
a set of two-dimensional constitutive equations and used them to model various planar load- 
ing cases, including cyclic loading. Their results compared quite well with experimental re- 
sults. 
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Zinov'ev and Sarbaev 21 recently used the endochronic theory to model composite laminate 
behavior. They used classical lamination theory in combination with endochronic equations 
to develop a material nonlinear laminate analysis code. They found good correlation with 
experiment for a variety of laminates subjected to in-plane loading. 


1.2.3 Finite Element Analysis with Material Nonlinearity 


Prior to the finite element method, stress analysis in the presence of inelastic material be- 
havior was a very intractable problem. Many approximate methods were developed to give 
engineers useful answers; however, the accuracy of these approximate methods left much to 
be desired. The advent of the finite element method and the digital computer provided engi- 
neers with the capability to obtain accurate solutions to complicated, nonlinear material 
problems. 

Inelastic numerical analysis has its roots in the simultaneous development of both the initial 
strain method 24 and the direct stiffness method 23 . In the initial strain method, the elastic 
equations of equilibrium are modified to include the inelastic strains which cause no change 
in stress. The marriage of the initial strain method to the finite element method was initiated 
by Gallagher and his Co-workers 25 2B . They used the initial strains to calculate an initial force 
vector (sometimes referred to as a psuedo force vector) which was added to the finite element 
force vector. 

In the direct stiffness method, sometimes called the tangent modulus method, the nonlinear 
problem is divided into a series of piecewise linear problems by stepping through the solution 
and modifying the stiffness at each step. The finite element tangent modulus equations were 
developed by Pope 27 , Swedlow and Yang 28 , and Marcal and King a . 
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Most present day inelastic solution techniques are based on one of these two methods. 
However, the efficiency, accuracy, and complexity of the present day techniques have all in- 
creased substantially. Inelastic analysis, geometric nonlinear analysis, and the combination 
of the two are now routinely performed using these present day techniques. 

1. 2.3.1 Endochronic Theory Applications 


The first application of the endochronic theory of plasticity in a finite element analysis was that 
of Valanis and Fan 30 . They used the initial strain method formulation in parallel with an 
incremental/iterative solution technique to calculate the stress and elastoplastic strain fields 
in a copper plate containing two symmetrically disposed edge notches. The plate was loaded 
cyclically in its own plane along its outer edges. Theoretical results were presented, but ex- 
perimental comparison was not provided. 

Watanabe and Atluri 31 used a different form of the endochronic theory, completely analogous 
to classical plasticity theory, in a finite element analysis to compare and contrast the results 
of Valanis and Fan 30 . The structure of the new endochronic equations allowed the use of the 
tangent stiffness method as opposed to the initial strain method used by Valanis and Fan. 
Again, theoretical results were not compared with experiment. 

1.2.3. 2 Composite Material Applications 

Applications of the finite element method to the analysis of composite material structures ex- 
hibiting inelastic behavior have been quite sparce. Adams and Miller *, Adams 33 , and Foye 
34 used the finite element method to analyze unidirectional fiber-reinforced composites mod- 
elled on the microscale. Their models envisioned arrays of elastic fibers embedded in an 
elastic-plastic matrix material. Such microscopic analyses can provide an understanding of 
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the inelastic behavior of unidirectional laminae; however, it cannot be extended efficiently to 
multidirectional laminates. 

Renieri and Herakovich 35 and Nagarkar and Herakovich 36 used a Ramberg-Osgood 
constitutive model in a quasi-three dimensional finite element analysis of interlaminar 
stresses in finite width laminates. Cross-ply, angle-ply, and two quasi-isotropic 
graphite/epoxy laminates were studied. Griffin, Kamat, and Herakovich 37 used a Hill-type 
yield criterion and an incremental plasticity approach in a fully three-dimensional finite ele- 
ment analysis of notched composite laminates. Material nonlinear isotropic hardening was 
modelled using the Ramberg-Osgood model and the constitutive equations were assumed 
temperature dependent. The incremental solution was performed with nonlinearities intro- 
duced as psuedo loads computed from initial strains. 

Using a previously developed elastic-plastic constitutive theory for composites 38 - 40 , Bahei- 
El-Din, Dvorak, and Utka 41 performed a finite element analysis of Boron/Aluminum plates of 
various layups with and without cutouts. Their analysis shows good agreement with exper- 
iment for uniaxial loading. Several other researchers also have recently applied various 
constitutive theories to inelastic finite element analysis of fibrous composites 42 - **. 
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2.0 The Endochronic Theory 


2.1 Introduction 


As noted in the literature review, macro-level nonlinear constitutive theories modeling the 
response of fibrous composite materials are limited in number and validity. The few theories 
available were considered, and due to its firm thermodynamic foundation and demonstrated 
correlation with experimental results, the endochronic theory 20 was chosen for this study. 
Pindera and Herakovich 20 have shown their anisotropic endochronic theory to model lamina 
response very well for in-plane loading. This study extends the anisotropic equations to three 
dimensions in order to analyze laminates with three-dimensional stress states. 

Vaianis 11 12 originally developed the endochronic theory to explain certain responses exhib- 
ited by ductile metals in the nonlinear range. Effects such as cross-hardening in tension due 
to torsion as well as the formation of hysteresis loops in loading/unloading cycles, not easily 
or accurately explained by classical plasticity theories, were accurately predicted by Vaianis 
for certain metals. The theory, based on irreversible thermodynamics, was developed using 
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the concept of internal variables. As previously mentioned, the theory is centered around a 
deformation scale which is assumed to be a material property. 

The deformation scale (also referred to as the time scale) can be either time dependent or 
time independent. In this study, the deformation scale is assumed to be time independent. 
Also, the deformation scale may be a function of stress or strain. This study uses a stress 
based scale (ala Pindera and Herakovich) to make the determination of endochronic constants 
more direct. This will become more evident later in this chapter. A primary feature of Pindera 
and Herakovich's development is the use of a multiple set of internal variables, which reduce 
the strain functions to power law approximations. Whereas metals are linear up to a yield 
point, facilitating a single internal variable approximation, epoxy resin composites generally 
show a gradual softening, with no well defined yield point. Multi-internal variables provide 
such behavior. 

Whereas Pindera and Herakovich used a single time scale z which is a function of all the 
stress components, Zinov'ev and Sarbaev 21 used multiple, separate time scales z, (i = 1 ,2,3) 
which were functions of a 2 , and t 12 , respectively. Valanis 18 similarly used the hypothesis 
of the existence of two types of intrinsic time for isotropic materials: the "hydrostatic" and the 
"deviator" types. This made it possible for Valanis to distinguish two deformation mechanisms 
in the inelastic range: the volumetric and the shear mechanisms, each characterized by its 
own intrinsic time scale. The use of multiple time scales allows for a more detailed material 
model and can, theoretically, yield more accurate analysis, at the cost of added complexity. 
The single time scale used by Pindera and Herakovich yielded good correlation with exper- 
iment for laminae, thus the added complexity of multiple time scales was felt to be unneces- 
sary for the present study. 
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2.2 Development of 3-D Endochronic Equations 


The endochronic equations are based on two fundamental laws of thermodynamics governing 
an irreversible process: 1) Conservation of Energy, and 2) The Dissipation (Clausius-Duhem) 
Inequality. From this base, Pindera and Herakovich 18 rigorously developed the equations for 
strain in terms of stress, a deformation scale z, and materia! constants. 


In the equations that follow, all stresses and strains are expressed in reduced notation as is 
standard in the composites literature, 


- 


" 


- 


- 

El 


e 11 




<*11 

s 2 


s 22 


o 2 


CT 22 

E 3 

— 

e 33 

and 

<*3 



a 33 

^4 


Y23 


<*4 


t 23 

e 5 


Yl3 


<*5 


t 13 

e 6 


Yl2 


<*6 


t 12 


( 2 . 1 ) 


where the 1-2-3 coordinate system corresponds to the material principal coordinates with the 
1-direction corresponding to the fiber direction (Figure 1). 

The reader is referred to Pindera and Herakovich's work 19 20 for a detailed thermodynamic 
derivation of the equations that follow. In the following equation set (2.2), the equations for 
normal strain (ej, i = 1,2,3) are obtained from a small modification of those developed in 
equation set (2.30) of Pindera and Herakovich's work 19 . The equations for shear strain (e„ 
i= 4,5,6) are identical to those shown in equation set (2.16) of Pindera and Herakovich. 

Pindera and Herakovich appropriately modelled the fiber direction response as a stiffening 
response because a, was always positive in their experiments. Because the laminates ana- 
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lyzed in this study are loaded in compression, the normal stress (a,) in the fiber direction is 
generally compressive. Therefore, rather than a stiffening response in the fiber direction, this 
study uses the softening response which is exhibited by compression o, - r,, experiments. This 
softening behavior is assumed to be uncoupled from the other stress components and is ac- 
counted for in the expression for e, by a new function f(o 1 ). The form of this function is defined 
later in this chapter. Notice also that the expressions for the nonlinear portions of e 2 and e 3 
(the integral quantities) in (2.2) are independent of a,. Equation set (2.2) is the complete set 
of three-dimensional equations needed for this study. 

As shown by Pindera and Herakovich, the anisotropic endochronic equations can be ex- 
pressed in the following form, where the term fjcr,) has been added to account for the softening 
as explained previously. For a single internal variable q, (see Pindera and Herakovich 19 for 
details on the internal variables qj, 


El = -^T + "«*i> 

— 5Gq r 2— 7 

E 2 = -fo 2 + Jo b 2 2 (z - z')c 2 (z')dz' + J 0 B 23 (Z - z')o 3 (z')dz' 

— <3Gq r 2 — 7 

e 3 = -35— + Jo b 22(z - z')c 3 (z')dz' + J 0 B 23 (Z - z')a 2 (z')dz' 


= 


-3G , 

d<iA 


0 + J 0 z E.a 4 (z')e-^ (z - z V 


65 " 15 ~ + Jo F * ^(z')e X,(Z ~ Z V 


e ® = ilr + J° F * CT6(z ' )e Mz-z ’ d2 ' 


6 


( 2 . 2 ) 


where 
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B 22 (z - z') - E-e-^ z “ z,) + 


B 23 {z-z')= -E*e ^ (z Z) + -^-*e X * (Z Z) 


E _ (B 22 B 23 ) 

2(b 22 — b 23 ) 


F = 





D — |-{[2B 12 b 11 — 2B 1 -|b 12 ] — P^(b 22 + b 23 ) 2B 12 b 12 ]a 1 } • (B 12 a 2 + -^-(B 22 + B 23 )] 


(1 — a 1 a 2 )[b 11 (b 22 + b 23 ) 2b 9 2 ] 



z' = Integration Variable 
L), X 2 = Eigenvalues of equation (2.23) in 19 

19 

o 1p a 2 = Normalized eigenvalues, equation (2.26) in 


In the above, dz = ^Si’dc^doJ’ where the terms of the fourth order, positive-definite tensor S,, 
are material parameters. Likewise, b^ , B tj , and Cy are all positive-definite fourth order tensors 
whose terms are material parameters. The X, terms are referred to as hardening exponents 
which reflect different dissipation modes. G 0 is a potential function defined as a stable equi- 
librium state, in directly measurable quantities (i.e. no internal variables). Extending 
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equations (2.2) to multiple internal variables requires a new definition for the B (J (z) terms in 
(2.2) and a new form of the expressions for e 4 , e 5 , and e 6 . 

The new expressions for e 4l e 5 , and e e with multiple internal variables qf (a = 1 to M) 19 are: 

c5Gq r z — 

e 4 = + Jo B 44<Z - Z')<T 4 (z')dz' 

“ 3Ga Z — 

C 5 + Jo B 66(Z - Z')CT 5 (2')dz' (2.3) 

5Gn 2 — 

e 6 = ~g~ + J 0 B 66 (Z - z')a 6 (z')dz' 


where all B^z) in (2.2) and (2.3) are now given as: 


B 22 (z) 


= ± y (B“ 2 + B23) 2 c -X?(z) , (B 22 -B“ 3 ) 2 -X?(z) 

2 a=l[ (b^ 2 - b^ 3 ) (b^ 2 -b^ 3 ) 


B2s(z) 


< B h + B 23) 2 C -X^(Z) (B 2 a 2 — B 2 3 ) 2 -x?(z) 

2 a = i|_ (b 2 2 — ^23) (b 22 -b^ 3 ) 


*a~1 (b 22 b 2 3) 


2 a - 1 b 66 


r a , r a 
a = ^22 ^ ^2: 

boo + bo< 




23 


23 


°22 

A-o — Ka — 


r a _ r a 
^22 u 2' 

Doo Do' 


23 


23 
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and 



M 


number of internal variables 


The expressions for e 2 and e 3 show the response in the direction transverse to the fibers to 
be governed by two independent sets of hardening exponents, A.f and A.?. Because A.J = A.§, 
coupling exists between transverse normal and shear response in the plane of isotropy. As- 
suming the hardening response is similar in shear and transverse tension in the plane of 
isotropy yields b? 3 = Bf 3 = Q 3 = 0. This simplifies the expressions for the B;j parameters. 



®22( z ) ~ 


“ (g^le-^lz) 

a— 1 bj2 


B 23 (z) = 0 


(2.4) 


® 44 < z ) ^^ 22 ^ z ) 


^66( z ) — 


l (^6 ) 2 c ~Xg(r) 

“ = 1 ^66 


This greatly simplifies equation set ( 2.2 ) to: 


®i = 


-gG 0 

<3(7., 


+ i{o,) 
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e 2 = 


E 3 " 


e 4 = 


e 5 ~ 


e 6 - 


-5G 0 
do 2 

+ Jo Z B 22 (z - 

z')cr 2 (z')dz' 

-3G 0 

+ Jo% 2 (z - 

z')a 3 (z')dz' 

-3G 0 

5ct 4 

+ yJo 2 B 22 (z - 

- z')a 4 (z')dz' 

-5g 0 

do 5 

+ Jo B 66( z “ 

z')a 5 (z')dz' 

~dG 0 

do 6 

+ Jo B 6 6 ( z - 

z')a 6 (z')dz' 


( 2 . 5 ) 


To obtain a power-law approximation for equations ( 2.5 ) , it is necessary to integrate the 
equations by parts. This is done for the expression e« below: 


e 6 “ 


~ dG 0 

do K 


+ Jo X — nf— -e x?(z Z ) <r 6 (z')dz' 
J ~ = i b 66 


e 6 " 



( 2 . 6 ) 


where 


U = c 6 (z') 


M (Brr) 2 _ 

dV = l — e ** (z z >dz' 
o=i b KR 


dU = 

dz' 
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Recalling that, 


J 0 z UdV = Uv|* - J z VdU 


we have 


~dG 0 

do* 
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X - Jo X-JJ 

a - 1 Ag bgg 0 a - 1 Ag 
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bee dz ' 
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However, 



^-[a 6 (z) - a 6 (0)] 


f z M 1 (Bg 6 ) 2 &r e (z’) 
■'°a = iXe bg 6 dz' 


therefore, 


,a .2 


— 3G n M 1 (Brr) £ 

c « - 1^- + <w» I tV .? « - • v ’> + 

0<J 6 01=1 Ae bee 


s x^d - ifv 
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Equation (2.7) corresponds to equation (5.17) in Pindera and Herakovich 19 . 
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Similar expressions hold for the remaining principle strains. Assuming a large number of 
internal variables (M) represent the material response, a power law approximation may be 
introduced as follows: 


1 


(B ? 6 ) 2 


b 66 


d 


- e-^ z ) s 


qO 

b 66 


z n * 


for n 6 > 1 


(2.8) 


where 


Beg is an arbitrary constant whose subscripts relate it to B 66 


n 6 is a constant whose subscript relates it to B 6 e 


The deformation scale z requires that a 6 (0) = cr s (0) = a 4 (0) = o 3 (0) = o 2 (0) = 0 . Using the 
power law type expression (2.8) for all the strain components gives the following form of the 
endochronic equations: 
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~dG 0 

3o 1 


+ fla,) 


-3G, 
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-+ jo B 22‘(2 ~ Z) ~^r dz 
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(2.9) 
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B °2 is an arbitrary constant whose subscripts relate it to B 2 2 
n 2 is a constant > 1 whose subscript relates it to B 2 2 


2.3 Incremental Form of Endochronic Equations 


The equations developed in the preceding section are not suitable for use in a finite element 
program. This is because the strains are functions of complex integrals involving the defor- 
mation scale z . Calculation of the parameter z requires knowledge of the load history. For 
simple cases, such as tensile or compressive loading of a lamina, the stresses are propor- 
tional to one another and these proportions do not change with increasing load. This allows 
z to be expressed as a proportionality constant times the desired component of stress. As 
an example, consider uniaxial tensile loading of a lamina transverse to the fibers. 

dz = N /S 22 da 2 dCT 2 = yj S 22 dcr 2 (2.10) 

z = \/ S 22 °2 

Similar expressions may be developed for other lamina orientations under a uniaxial load. In 
fact, the endochronic constants are calculated using these uniaxial load cases. This is shown 
later in section 2.4. 

In a laminate, the stresses need not be proportional to one another. Since the load history 
of the stress components is not known a priori, equations ( 2.8 ) may not be integrated exactly. 
However, equations (2.8) may be integrated incrementally by assuming the stress compo- 
nents to be proportional to one another over each load increment. This allows the stress 
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component load histories to be nonproportional in a piecewise sense. Consider the ex- 
pression for e«: 


e 6 - 


~ dG 0 


+ JX< Z - 2 ’ ,n ‘'^r dz ' 


Written incrementally, this becomes 


e 6 = + Jo 1B 66( z n " O^dz' + Jo’bJ^Zn - z')" 6 S6 2 dz' + 


+ Jq N ^66( z n z') ne S6 N dz' 


or 




(n 6 + -jj {S6 2 [( z n Zi) 9 (z N z 2 ) n ® ^ + + S6 n (z n 
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dz 


Aa 6 
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k = 1,2, N 


z N ~ 2 Az k 

k = 1 


Az k = ^/SjjAa^Aaf 


N = number of increments 


( 2 . 11 ) 
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The integration accuracy is a function of the increment size. The smaller the increment, the 
more accurate the integration. In the limit, as Az goes to zero, the integration becomes exact. 
For mildly non-proportional load histories, the error introduced by this process should not be 
great. 


2.4 Determination of Endochronic Parameters 


Determination of the endochronic constants requires analyzing many test results. Burns et 
al 1 performed all the laminate tests as well as the 10° and 45° off-axis tests and the uniaxial 
[0] and [90] compression tests. All losipescu tests as well as the 15° and 30° off-axis tests 
were performed by Jack Beuth 50 at VPI. Data from at least two and in some cases three tests 
were used in determining the constants for this study. The final value of a particular 
endochronic constant is the average of the values determined from all tests. Any exceptions 
to this procedure are noted. Appendix A shows the least squares method of representing a 
"best fit" line through the experimental data. 

The endochronic equations may be written in terms of reversible strain (e?) and irreversible 
strain (e| r ). The reversible strains are represented by derivatives of the potential function G 0 
in equations (2.9), whereas the irreversible strains are represented by the integral quantities. 

e, = e, R + e, IR (2.12) 


where 


R ~dG 0 
1 fo, 


i = 1.2. 
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and 


The Endochronic Theory 


22 
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e!| R = f<<?i) 

+ 

p ,R - r 

E 2 Jc 

;B° 2 • (z - z') n ’-^rdz' 


JR = r : 

e 3 Jc 

;B° 2 -(z-z') n *^dz' 

• 

CO 

II 

^•(z-zy^ltdz 

• 

JR _ f : 

E 5 Jc 

;B° e .(z-z') n *-|fdz' 


JR = r 

E s J c 

;B° e -(z-z') n# -^fdz' 


2.4.1 Reversible Response 

The potential function G 0 has not been discussed in any detail to this point because its form 
is strongly material dependent and stress sign dependent. As noted earlier, Pindera and 
• Herakovich used a function G 0 which modeled stiffening behavior in the fiber direction of 

C6000/PMR15 graphite-epoxy under uniaxial tension. AS4/3502 exhibits similar tensile be- 
havior. However, because the laminates in this study are loaded in compression, the fibers 
^ are in compression except, possibly, in regions near the hole. A typical uniaxial compression 

test of AS4/3502 1 exhibits softening as shown in Figure 2. This softening behavior is modeled 
in this study as an uncoupled event. Therefore, it is not represented by a function of z nor is 
it included in the potential function G 0 , ala Pindera and Herakovich. Instead, it is expressed 
as a separate function f(a,) which is determined in the following section. 
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The function G 0 needs to model only reversible linear response and is defined as 


-G 


0 ” _ 2‘ A ij a l a J 


(2.13) 


Therefore, 


R dG 0 

1 d<J t 


A„a, 


i = 1,2, 3, 4, 5, 6 


(2.14) 


A thermodynamic constraint on G 0 requires the fourth order tensor A M to be positive 
definite 19 . 


|im g2 ( -G o) - A^Sa, ^ 0 (2.15) 

CTj -»CO 

The fourth order tensor Ay is the compliance matrix. It represents the initial compliance of the 
material and is determined from appropriate tests. 


2.4.2 Response Along the Fiber Direction 

The function f(o,) originally introduced in equation (2.2) is determined from a uniaxial 
compressive test (Figure 2) with the fibers aligned along the specimen's axis. A simple power 
law fits the experimental data very well; therefore, the functional form of ffo,) was chosen to 
be; 


f(a 1 ) = ± B.,., la, I" 1 where n 1 > 1 (2.16) 

The absolute value sign is used because Oi is usually negative. The function f(a,) is negative 
if o, is negative and positive if a, is positive. 
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cr 1 (KSI) 
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Substituting equations (2.13) and (2.16) into equation set (2.8), the complete form of the 
endochronic equations can be written as 


£^ — A.jjOj i 8^ |<J 1 1 




=2 - *2|°| + Jo B °2 ' < Z - -dz' 


^3 


. _ , fZ«0 . _,*n, Sa 3 {z') 
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e 4 = + |j 0 Z B^-(z - zf^dz' 


(2.17) 


e 5 = AsjOj + JoBge • (z - z') n « - ^ 1 dz' 


e 6 


= Aejcrj + J 0 Z B° 6 -(z - z')" 6 — ^ 7 —dz' 


2.4.3 The Compliance Matrix 


The reversible strain of the material is expressed in terms of stress and the compliance matrix 
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(2.18) 
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Under the assumption of transverse isotropy in the 2-3 plane, A 33 = A 22 , A 13 = A 12 , and 
A 55 = A ge . This simplifies the above expression to 
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(2.19) 


Equation (2.19) contains only five independent constants (A,,. A 12 , A K , A u , and A se ). The addi- 
tional A 23 constant is a function of A 22 and A M (i.e. 2(A 22 - A 23 ) = A u ). 


2.4.3.1 The A 11 Constant 


The constant A„ is determined by considering the compressive loading case of a, * 0 and 
all CT| = 0, (i — 2, 3, 4, 5, 6 ). The reversible strain in the fiber direction is then 

e? = A^o, (2.20) 


and 



A n 


However, from (2.17) 


5ei dz] j 

I a, =0 
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Therefore, 


An ~ 


dz] 

5<Ti 'ct, =0 


( 2 . 21 ) 


The data from two compression tests on unidirectional [0) AS4/3502 were essentially linear in 
the range of o, = 0 to a, = 25 ksi. A least squares fit through these points gave a straight line 
with a slope m (or E^. An is the inverse of m (or Ei). 

Averaging the values for A„ from both tests (Table B.1(a)) yields: 

A^ = 5.12687 x 10“ 8 (psi -1 ) (2.22) 


2A.3.2 The A n Constant 

Like A 11t the endochronic constant A 12 is also calculated from a [ 0 ] uniaxial compression test. 
Consider the expression for as shown below: 

e 2 = A 12 Ci for a., * 0 Oj = 0 (i = 2,3, 4,5, 6) (2.23) 


Now 


dz* _ do i 
3e 1 12 dzi 


and 

8 ±\ _ A | 
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OZ y '<y, =0 ^ 'a, =0 


From (2.21) 
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and by definition 


fog 

dEj 


< 7 1 =0 


= v 12 


Therefore, 


A12 — A^v 12 


(2.24) 


In order to calculate A 12 , v 12 must first be determined. This is done by fitting a least squares 
line through a plot of e 2 versus e, (Figure 3). The least squares fit is used in the range of data 
points from e 1 = 0 to e 1 = 0.13%. The value of v 12 is the slope of this line. 

Averaging the values of v 12 from both tests (Table B.1(b)) yields 


v 12 = 0.333 


(2.25) 


Combining (2.22), (2.24), and (2.25) gives 

A 12 = -1.70722 x 10 -8 (psi -1 ) (2.26) 

2.4. 3.3 The A 22 Constant 

A 90° uniaxial compressive test is used to determine A 22 . In this test, the fibers are oriented 
transverse to the loading axis. Because all o, (i = 1,3, 4, 5, 6) are zero, the expression for e£ is: 

e 2 = A 22 a 2 (2.27) 
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and 


de* <?e 2 | 

3a 2 22 3 ct 2 I O2=0 

Utilizing the plots of o* versus e 2 shown in Figure 4, a least squares fit is used in the range 
of data points between o 2 = 0 and a 2 a 8 ksi. The endochronic constant A a is the inverse of 
the slope of this line (modulus E 2 ). 

Averaging the value of A M from both test runs (Table B.1(c)) yields 

A 22 = 6.55133 x 10 7 (psi _1 ) (2.28) 

2A.3.4 The A ^ Constant 

The endochronic constant A44 is determined from an losipescu shear test in the 2-3 plane. This 
type of test induces an approximately pure 2-3 shear stress o 4 (t 23 ) in the central region of the 
specimen 47 . Figure 5 shows the test specimen and idealized loading conditions used for this 
test. (The symbol 0 in this figure is only valid for 1-2 plane tests. For 2-3 plane tests, the fibers 
run in the direction of page depth). Figure 6 shows plots of a 4 (x 23 ) versus e 4 (y 23 ) for two tests. 
The endochronic constant is determined by fitting a least squares line through the initial por- 
tion of the curve for the range of data points between o 4 = 0 and o 4 = 6 ksi. The constant 
A* is the inverse of the slope of this line (modulus G*). 

Averaging the values of A* from tests 1 and 2 (Table B.1(d)) yields 

A44 = 2.18049 x 10~ 6 (psi -1 ) (2.29) 
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Figure 4. o 2 Versus c 2 for 90° Uniaxial Compression Tests of AS4/3502 
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FJ ,f| = Distributed Forces 
p s Fiber Orientation Angle 


Figure S. Test Specimen and Loading Conditions Used for losipescu Tests 
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Figure 6. cr 4 Versus e 4 for losipescu Tests In the 2-3 Plane of AS4/3502 
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2.4.3. 5 The A 2 3 Constant 


Because a transversely isotropic material contains only five independent constants, A 23 may 
be determined indirectly. Consider the following expression for the shear modulus G y 
(i , j = 1,2,3) in a plane of isotropic symmetry. 


G„ = - 


2(v,j + 1) 


(2.30) 


G 23 = 


E 2 

2(v 23 + 1) 



2 < V 23 + H = 2 / v 23 , J_\ 
E 2 \ E 2 E 2 ; 


However, by definition, 


A44 - 




A 23 “ 


v 23 

E 2 


Therefore, as shown earlier, 


A 44 - 2(A 22 — A 23 ) 



From the previously determined values of A K and A M , A 2J is calculated to be 
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2.4.3.6 The 4 66 Constant 

The endochronic constant A 6 , is determined in an analogous manner to A*. An losipescu 
shear test in the 1-2 plane yields a region of approximately pure a, (t 12 ), from which a 6 versus 
e s (Yi 2 ) way be plotted. Figure 7 shows such a plot for two tests. A least squares fit line 
through the initial portion of each curve yields a slope m (modulus G 12 ), the inverse of which 
is Agg. The equation for e£ under the condition o t ^ 0, a, = 0 (i = 1,2, 3,4,5) is shown below. 

e 6 = ^66^6 (2.32) 


and thus 

^6 . dE6 | 

do 6 66 du 6 l 0g = 0 

Because the test data was available, off-axis tests of 10°, 15°, 30°, and 45° were also utilized 
to determine Ag 6 . Figure 8 shows a 8 versus c, for the four off-axis tests. 

Averaging the values of A 66 from the above six tests (Table B.2(a)) yields 

A 66 = 1.18373 x 10' 6 (psi -1 ) (2.33) 

2.4.4 Irreversible Endochronic Constants 


The irreversible endochronic constants are determined by plotting In e[ R versus In a, for vari- 
ous lamina tests. These plots should be a straight line because when only one stress com- 
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ponent (a,) is present, the irreversible portion of strain (eJ R ) is represented by a power law 
(2.12). Symbolically, we can represent the power law form of e! R as follows: 


Ip |> 

e ( = CjCT|' (i = 1,2, 3, 4, 5, or 6, no sum on i) 

where 

C ( = functions of the irreversible endochronic constants 

ki m ni 
k 2 = n 2 + 1 

k 3 = n 2 + 1 

k 4 = n 2 + 1 

k 5 = n 6 + 1 

kg = n 6 + 1 

Taking the logarithm of both sides of the above equation gives: 

In e| R = kj In Oj + In C, 

In the classical x-y coordinate system, the above can be represented as: 

y = mx + b 

where 

y s in e| R 
x = In Oj 


(2.34) 


(2.35) 
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In this coordinate system, m is the slope of the line and b is the y-intercept. The above rep- 
resentation is used in evaluating all of the irreversible endochronic constants. 

2.4.4.1 The n and B Parameters 

The same tests used to determine the Ai, parameters are utilized in calculating the 
endochronic constants appearing in the irreversible strain expressions. These constants in- 
clude n, , B 11t n 2 , BJ 2 , n„, and Bg 6 (2.17). 

iii and Bn 

The constants n, and B„ are calculated from the 0° uniaxial compression tests. The ex- 
pression for e',* is written as follows: 

pIR = R rr n 
e 1 


Using equation (2.34), 


c i ~ B ii 
k i = n i 

For a power law fit, this plot should yield a straight line with slope k, and y-intercept In C,. 
Figure 9 shows the experimental curves to be best represented by a straight line in the range 
of data points from In o, = 11 to In o, = 11.7. 
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Figure 9. In e 1 ir Versus In a< for 0° Uniaxial Compression Tests of AS4/3502 
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Averaging the results from tests 1 and 2 (Table B.2(b)) yields 


n, = 2.44791 


(2.36) 


Bn = 1.90645 x 10 


(2.37) 


n 2 and B ° 22 

The constants n 2 and B£> are determined from the 90 0 uniaxial compression tests. The ex- 
pression for 4 r (2.17) is 


where 


JR 

e 2 


= 


22 


<Z - zrf&PU? 

dz' 


dz = ^SjjdOidOj 


(2.38) 


Under pure uniaxial loading, a 2 * 0 and o, = 0 (i = 1,3, 4, 5,6). The expression for dz simplifies 
to 


Therefore, 


dz = ■ v /s 2 2dtf2<fo2 = V S 22 d<T 2 




22 a 2 


8a 2 _ -) 


(2.39) 


and 


The Endochronic Theory 


42 


Integrating, 


JR 

e 2 


= J 0 B22-(2-Z') n * 


1 

V S 22 


dz' 


pO 

JR _ B 22 

E 2 /==■ 

22 


V & 


(z - Z') n * +1 


n 2 + 1 


pO n 2 + 1 
B 22 z 


V S 22 < n 2 + 1) 


JR _ ®22 (®22) 2 n, +1 
-Co* 


e 2 = 


n 2 + 1 


Using equation (2.34), the above can be symbolically represented as: 


JR 


s C 2 a 2 * 


where 


C 2 


B 


22 " 


nj 

(S 22 ) 2 
n 2 + 1 


1^ — n 2 + 1 


(2.40) 


In the above equations, In C 2 is the y-intercept and k* is the slope of the In e^ r versus In o 2 
curve. 


The plots of In versus In a 2 for two tests are shown in Figure 10. A straight line is best re- 
presented over the range of In a 2 = 9.9 to In <j 2 = 10.3. 

Averaging the values of n 2 from both tests (Table B.2(c)) yields 
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In o 2 


Figure 10. In e 2 IR Versus In a, for 90° Uniaxial Compression Tests of AS4/3502 
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n 2 = 2.45967 


(2.41) 


It was necessary to calculate n 2 before B 22 because is a function of n 2 as well as S 22 (2.40). 
Because S a has not yet been determined, Bg 2 will be left for now as a function of S„. 


Since In C 2 is the y-intercept, we have 


In 


( 


ri_ 

622 (^ 22 ) 2 \ 
n 2+ 1 ) 


y intercept 


2 = e y ln t.rc.pt 

n 2 + 1 
0 ““ 

B 22 = (S 22 ) 2 (n 2 + 1 ) e 71 "**™* 

Averaging the values of Bg 2 from both tests (Table B.2(d)) yields 

B 22 = s 22 1 ' 22984 (2. 3700 96 x «T 18 ) 


(2.42) 


n« and B 0 « e 

The constant n 8 is evaluated in a similar manner to n 2 using the two losipescu tests in the 1-2 
plane along with the 10°, 15°, 30°, and 45° off-axis tests. The expression for ei* is written as 
follows: 


JR _ D 66V°66J * _n, +1 
£s r)«” i 06 


Using equation (2.34), the above can be symbolically represented as 


JR /s kt 

^6^6 


(2.43) 
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where 


o 

c s B 66 (S 66 ) 2 
6 n 6 + 1 

ke = n 6 + 1 

The irreversible strain ei R can be written as: 

e 6 R = e 6 — e 6 = e 6 — ^66 a 6 

The irreversible strain e|, R may be calculated in two different ways. The first scheme uses the 
average value of A 66 determined in the previous section. The second scheme uses the value 
of A w determined for each separate test. Since it is not known a priori which scheme will yield 
the most accurate value of n 6 , both schemes are used to determine separate values of n 6 . 
These values are compared with experimental results later to determine a best fit value. 
Figure 11 shows a least square fit to be best represented in the following ranges: losipescu 
Test 1 (In o 8 = 8.7 to 9.5), losipescu Test 2 (In cj 8 = 8.62 to 9.51), 10° Off-Axis Test (In a, = 8.39 to 
9.16), 15° Off-Axis Test (In a 6 S8.14 to 9.0), 30° Off-Axis Test (lno e S8.62 to 8.91), and 45° Off- 
Axis Test (lno 8 = 8.25 to 8.46). The following result is based on the average A 88 value of 
1.18373 x 10" 8 (ps'r 1 ). 

Averaging the value of n 8 from the six shear tests (Table B.3(a)) yields 


n 6A = 2.54504 (2.44) 

where the subscript 6A represents the value of n 8 calculated using the average value of A 86 . 
Averaging the value of n 8 from the six shear tests (Table B.3(b)) yields 


n 6B — 2.53355 


(2.45) 
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Figure 11. In e, IR Versus In CT, for losipescu Tests and Off-Axis Tests of AS4/3502 
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where the subscript 6B represents the value of n s calculated using the individual test values 
of A w . 

Comparing the plots of In ej, R versus In a 6 in Figure 11 shows that the 0 0 losipescu test best 
represents a straight line. Therefore, the value of n 6 from this test is also considered in the 
best fit scheme in the following section. Actually, two values of n t exist for the 0 0 losipescu 
test; a value based on the average A M compliance and a value based on the A S6 compliance 
of the 0° losipescu test. These values are represented by the subscripts 6C and 6D respec- 
tively and are listed below. 


n 6C = 3.04940 

(2.46) 

ngp = 3.1854 

(2.47) 


Because both Bg 6 and S M must be calculated from a pure shear test in the 1-2 plane, one of 
them is redundant (see equation (2.43)). Since Bg e is simply a multiplicative constant which 
shifts the curve of In eJF versus In a, to the left or right, it is arbitrarily set equal to the value 
of 1 in this work. 

Bee = 1 (2-48) 


2.4.4.Z The Sy Parameters 

With the exception of S* and Sgg, the S i( parameters are calculated from the off-axis coupon 
tests. S44 is calculated from losipescu tests in the 2-3 plane and is calculated from losipescu 
tests in the 1-2 plane. 
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The Constant S 44 


The equation governing the irreversible strain response in a 2-3 plane losipescu test is shown 
below. 




(2.49) 


In the same manner as equations (2.38), (2.39), and (2.34), the above can be symbolically re- 
presented as 

jr _ ~ k 4 

e 4 “ C 4 ct 4 


where 


C A = 


cjO 

°2Z 


n 2 

044 )T 
2(n 2 + 1) 


■<# 


+ 1 


k 4 = n, + 1 


Because n 2 has already been calculated and Bg 2 has been calculated as a function of S 22 , S M 
may be calculated in terms of B| 2 . However, examining the stress-strain behavior of the 2-3 
plane losipescu tests as shown in Figure 6, the stress-strain response is essentially linear \ 
This observation means the irreversible portion of e 4 is zero. If ei R is zero and n 2 and B§ 2 have 
been determined (and are used in the irreversible strain expressions of other strain compo- 
nents), then S*, must be set equal to zero. 


S44 = 0 


(2.50) 


1 These tests did not fail in a shear mode. Therefore, their ultimate strain values may be significantly 
larger. It is possible that these specimens would show non-linear behavior at higher strain levels. 
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The Constant S 23 


In the same manner as the A,, parameters, the S (j parameters are assumed to obey 
transversely isotropic conditions. Therefore, S 4j (j = 1,2, 3, 5,6), S 5j (j = 1,2, 3, 4, 6), and S 8j 
(j = 1.2, 3,4, 5) are all set to zero. Likewise, S 13 = S 12 , S 33 = S^, = S 66 , and in general, 

Sjj = Sjj. The expression for S 23 in terms of S 22 and must hold as well. 


S23 - S22 



(2.51) 


However, because S u = 0, S 23 must equal S 22 . 

S 23 = S 22 


(2.52) 


The Constant S 4 « 

As mentioned previously, the S 6e parameter is determined from losipescu tests in the 1-2 
plane. The equation governing the irreversible strain state ej, R in this test was shown previ- 
ously to be 

JR _ ~ k 6 
e 6 — CgOg 8 


where 

0 -2s. 

= B 66 (S 66 ) 2 

Ce - n 6 — 
kg = n 6 + 1 

Because B§, = 1, the expression for S M may be written as 
The Endochronic Theory 


In 


/ (s 66 rr \ = 

\ n 6+ i ) 


yintercept 


(See) 2 = (n 6 + 1)e yinWre “>' 


S 66 = [( n 6 + 1)e y '"“"*P‘J n, 


Note that S M is a function of n 6 . Therefore, there will be four values of S 66 corresponding to the 
four values of n 6 (n^ng^n^, n 6D ). The following value of Sgg is calculated using n w . 

Averaging the values of S 66 from tests 1 and 2 (Table B.3(c)) yields a value for S e6A of 

S 66A = 1-96639 x 10" 13 (2.53) 

In an analogous manner, S MB , S^, and S 6SD are calculated to be 


Sgea = 1.805796 x 10 -13 

(2.54) 

S 66C = 1.26342 X 10" 12 

(2.55) 

S 66D = 1.839071 x io“ 12 

(2.56) 


The remaining S (j parameters are calculated from off-axis tests. 

The Constants S,,, S 12 , and S 22 

As mentioned previously, an off-axis test introduces not only the desired shear stress, but 
normal principle stresses parallel and perpendicular to the fibers. Therefore, an off-axis 
coupon test will have non-zero stress components c^, a*, and a, (t 12 ). The ratio of these 
stresses depends on the off-axis angle, but throughout the loading, equilibrium requires that 
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they be proportional. In other words, the ratios of these stress components do not change 
during loading (i.e. — = const., = const., etc.). This is a very important result because 
it allows the three remaining 8„ parameters, S 11t S 22 , and S 12 to be calculated from three 
separate off-axis tests *. 


The expression for the irreversible strain in the 1-2 plane (e}p) is written 


r IR -f 2 pO./, &t 6 <z') ^ , 
e e “ Jo B 66 • (z - z ) ~ dz , dz 


where for the combined stress state a„ a t , cr 8 

— V S ii dCT i S 2 2dcT2 + 2S 12 da 1 da 2 + Sggdog 

Writing do, and do 2 in terms of do 8 gives 

do, = cot <p da 6 
do 2 = tan <p da 6 

where 


(2.57) 


<P = Angle of Fibers from the Loading Axis 


Therefore, 


V®ti co * ‘P + 2S 12 + Spg + S 2 2tan 2 <p dag 
2 " V s ii cot + 2S 12 + Sg 6 + S 22 tan 2 <p a 6 


rtXsTto S,» rl S * nd ° ch,0 " l< ’ *•« written for strain in ton.. 

conmonU? ,S “re SSSw S « 1 para "' e,ars N °'« that «» °«-axis tens, the strain 
ticuiuo * would b * ^ * 
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Substituting the above expression for z into the integral equation (2.57) gives 


_ u 

IR bge r a Pg 

e 6 = ' n ‘ + 1 [ s ii cot <P + 2S i 2 + + S 22 tan 2 (p]Toe* +1 


(2.59) 


Using equation (2.34), equation (2.59) can be symbolically represented as 


JR _ ~ k. 
e 6 ^ 6 ^ 6 ° 


where 


C« ■ 


b 66 
n R + 


J [ S * 11 cot 2 (p + 2(S 12 + S 6 e) + S 22 tan 2 (p] 2 


ke - n 6 + 1 * 

For three off-axis tests of different angles, the above expression yields three equations for the 
three unknown parameters S 11t S a , and S 12 . To use these equations, the In ejp is plotted 
against the In a t for each off-axis angle. Figure 11 shows these plots for the angles of 10°, 15°, 

30°, and 45° (because four tests were available, four sets of three different off-axis tests were 
utilized). 

Pindera and Herakovich " rearranged the above equations to yield simple equations for de- 
termining 8,,. S M , and S 12 . These equations (eqs. (5.22) and (5.23) in Pindera and Herakovich) 
are shown below. 

2 

S 11 = (n 6 + 1)n 8 x Q 

S 22 = (n 6 + l)"^" x y 

S 12 = y {[( n e + ^Ase^i)] 7 ^ ~ [s^ctpi + S 22 t<Pi] - S 66 \ (2.60) 
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where 


Cl = 


[a66(<Pi) n ® 

* X' 
~ A 6e ((p 2 ) n « _ 

[t<Pi ~ t<p 3 ] 

- 

' * X 

_Ag6(<Pl) n * - 

* -i-1 

' A 66 ((p 3 ) n$ J 

[tip, - t<p 2 ] 


[c<p-, - c<p 2 ] [t<p 1 - tlp 3 ] 

- 

[c<p 1 - C(p 3 ] [tip, - t(p 2 ] 



} 


'-{■ 


r x jl 

LA 6 6(<Pi) n ® ~ 

• X' 

' A 66 (<p 2 ) _ 

[cip, - c<p 3 ] 

- 

' . X 

.Ag0(ip,) n 6 

* xl 

- A 66 (ip 3 ) n 6 J 

[cip, - cip 2 ] 


[tip, - tip 2 ] [cip, - C(p 3 ] 

- 

[tip, - t<p 3 ] 

[cip, - c<p 2 ] 


} 


* 

Aigg(<p) = — ; [Si-i cot 2 <p + S 22 tan 2 (p + (2S 12 + S 66 )]T 

n 6 + i *- 


where i can be either 1,2, or 3, tip, = tan*ipj and c<p s = cot 2 <p ( 

Using the above equation set along with plots similar to Figure 11 allows S„, S a , and S 12 to 
be calculated. Using the four off-axis tests available in conjunction with the four sets of n 6 and 
S 66 produces sixteen different sets of parameters. These sets are shown in Table 1. 


2.4.5 "Best Fitting" the Parameter Sets to Experiment 


It was mentioned earlier that a condition imposed on S if by thermodynamic constraints is it 
must be positive definite. Therefore, of the sixteen parameter sets shown in Table 1, only 
those which are positive definite are possible choices. Appendix C shows the conditions 
necessary for a tensor to be positive definite. From these conditions, Table 1 lists the pa- 
rameter sets which are positive definite and those which are not. The parameter sets which 
are positive definite are numbered in Table 2 for ease of reference. 

To determine which of the positive definite parameter sets best fits experimental test results, 
the parameter sets are used analytically in a comparison with off-axis test results and are then 
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Table 1. Possible S„ Parameter Sets Marked As Positive Definite or Non-Positive Definite 


Angle Set 
(15 °,30 °,45 8 ) 

Angle Set 
(10 °,30 °,45 °) 

Angle Set 
(10 °,15 °,45 °) 

Angle Set 
(10 °,15 °,30 °) 

2.54504E-00 
1.96639E-13 
S t1 2.34704E-14 

S 22 3.79959E-13 

S 12 -6.36714E-14 

njA 2.54504E-00 

1.96639E-13 
S„ 3.46034E-15 

Sa 3.19929E-13 

S„ -2.36515E-14 

n M 2.54504E-00 

S 66A 1.96639E-13 

S„ -8.53148E-15 

S K -6.57692E-14 

S 12 1.75194E-13 

n 6A 2.54504E-00 

S 66A 1.96639E-13 

S„ -1.07261E-14 

S 22 -1.04893E-12 

S 12 2.25771E-13 

Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

n 6B 2.53355E-00 

S 66B 1.80580E-13 

S„ 1.05681E-14 

S 22 2.63695E-13 

S 12 2.26620E-15 

n 6B 2.53355E-00 

S 66B 1.80580E-13 

S,, 9.79461E-16 

S 22 2.34929E-13 

S 12 2.14432E-14 

n 8B 2.53355E-00 

S MB 1.80580E-13 

S„ -4.76695E-15 

S 22 5.01046E-14 

S, 2 1.16729E-13 

n 6B 2.53355E-00 

S 6SB 1.80580E-13 

S„ -5.81861E-15 

S 22 -4.21020E-13 

S 12 4.40965E-13 

Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

ngc 3.04940E-00 

S 66C 1.26342E-12 

S„ 1.49223E-13 

S 22 2.86855E-12 

S 12 -4.38706E-13 

n 6C 3.04940E-00 

S 66C 1.26342E-12 

S„ 2.26566E-14 

S 22 2.48885E-12 

S 12 -1.85573E-13 

n 6C 3.04940E-00 

S 66C 1.26342E-12 

S,, -5.31935E-14 

S 22 4.92545E-14 

S 12 1.07215E-12 

n 6c 3.04940E-00 

S 66C 1.26342E-12 

S„ -6.70750E-14 

S 22 -6.16938E-12 

S 12 4.39207E-12 

Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

Non-Positive Definite 

n 6D 3.18540E-00 

S 66D 1.83907E-12 

S„ 1 . 1 0652E-1 3 

S 22 3.57840E-12 

S 12 -7.48855E-14 

n 60 3.18540E-00 

8 m 1.83907E-12 

S„ 1.17136E-14 

S 22 3.28159E-12 

S 12 1.22990E-13 

n 6D 3.18540E-00 

S 660 1.83907E-12 

S„ -4.75791E-14 

S 22 1.37453E-12 

S 12 1.10616E-12 

n 6D 3.18540E-00 

S 66D 1.83907E-12 

S,, -5.84304E-14 

S 22 -3.48663E-12 

S 12 1.35624E-12 

Positive Definite 

Positive Definite 

Non-Positive Definite 

Non-Positive Definite 
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Table 2. Positive Definite Parameter Sets 


Parameter Set 


1 


2 


3 


4 


5 


S ii 


1.49223E-13 


1.71360E-14 


2.34704E-14 


1.10652E-13 


1.05681E-14 


S 22 


2.86855E-12 


3.281 59E-1 2 


3.79959E-13 


3.57840E-12 


2.63695E-13 


S 12 


-4.38706E-13 


1.22990E-13 


-6.36714E-14 


-7.48855E-14 


2.26620E-15 


S 66 


1.26342E-13 


1.83907E-12 


1.96639E-13 


1.18391E-12 


1.80580E-13 


n « 


3.04940 


3.18540 


2.54504 


3.18540 


2.53355 
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used in a laminate analysis to compare with angle-ply laminates. The parameter set which 
best fits the experimental data is chosen. 

2.4.5.1 Off-Axis Test Comparison 

The expression for e 8 in an off-axis specimen (2.56) is 

B° n 6 

e 6 = A 66 a 6 + cot2( P + ( 2S i2 + s 66) + S 22 tan 2 (p]T'CTe 6 + 1 (2.61) 

From this expression, e 6 may be plotted against cx 6 for any desired off-axis angle (p. A simple 
program was written for equation (2.58) that uses S,,. S 22 , S 12 , S 68 , n 6 , and «p as input param- 
eters and produces data for <r 6 versus e 6 . These data were plotted against experimental data 
and the plots are shown in Figure 12 through Figure 16 for the various parameter sets. 

Parameter set 2 (Figure 13) fits the experimental data most accurately in an average sense. 
Parameter sets 1,3,4, and 5 all seem to match the 15°, 30°, and 45° experimental data very 
well. However, they do a poor job of modelling the 10° off-axis behavior. This is to be ex- 
pected because, as Table 1 shows, these parameter sets were all derived from the combina- 
tion of 15°, 30 c , and 45° off-axis tests. Parameter set 2, however, was derived from the 
combination of 10°, 30°, and 45° off-axis tests. Although it matches these tests well, parameter 
set 2 also does an adequate job of modelling the 15° off-axis test. 

2.4.5.2 Angle-Ply Laminate Comparison 

Unlike the off-axis lamina tests, the loading within the layers of a laminate is not proportional. 
Therefore, no simple expressions exist for strain in terms of stress. Instead, the endochronic 
equations must be integrated throughout the laminate's load history. This is done in a step- 
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Figure 12. Analytical/Experimental Comparison of Off-Axis Tests Using Parameter Set 1 
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(a) 10° Off-Axis Comparison 
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Figure 13. 


Analytical/Experimental Comparison of Off-Axis Tests Using Parameter Set 2 
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ct 6 (ksi) _ o;(ksi) 




^6 (%) ^6 (%) 


(c) 30° Off-Axis Comparison (d) 45° Off-Axis Comparison 
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Figure 14. Analytical/Experimental Comparison of Off-Axis Tests Using Parameter Set 3 
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0;(KSI) <X(KSI) 



(a) 10° Off-Axis Comparison (b) 15° Off-Axis Comparison 
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ct v (ksi) o;(ksi) 



€y(%) € y (%) 

(a) [( ± 10) 12 ], Laminate Comparison (b) [( ± 20)„], Laminate Comparison 



(c) [( ± 30),j] t Laminate Comparison (d) [( ± 45) 12 ], Laminate Comparison 
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Figure 17. Analytical/Experimental Comparison of Angle-Ply Laminate Tests Using Parameter Set 
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o;(ksi) cr v (Ksi) 




(c) [( ± 30),j], Laminate Comparison 
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Figure 18. Analytical/Experimental Comparison of Angle-Ply Laminate Tests Using Parameter Set 
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°y (KS,) a- y (KSi) 



€ y (%) € y (%) 

(a) [( * 10>«]. Laminate Comparison (b) [( ± 20)*]. Laminate Comparison 



(c) [( ± 30)„], Laminate Comparison (d ) [( ± 45)*]. Laminate Comparison 
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0;(KSI) <X V (KSI) 
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Figure 21. 


Analytical/Experimental Comparison 
5 


of Angle-Ply Laminate Tests Using Parameter Set 
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wise fashion and for each step, the loading is assumed to be proportional. The step-iterative 
constant stress and constant strain algorithms developed in Chapter 3 were used in a finite 
element program to analyze laminate problems using a four element mesh. The flow charts 
for this program are analogous to those in Figure 30 and Figure 32 in Chapter 3. 

The output from the finite element analysis was plotted against the experimental data 1 for 
C ± 103., [ ± 20]„ [ ± 30]„ and [ ± 45], angle-ply laminates. These plots are shown in 
Figure 17 through Figure 21. Examining these figures shows that parameter set 2 
(Figure 18) also fits the laminate experimental data most accurately. In fact, the correlation 
may be described as excellent. Many of the other parameter sets match one or two of the 
laminate tests well, but caused the solution algorithm to diverge for the other tests. This di- 
vergence is caused by certain stress components decreasing in magnitude (unloading) and 
is discussed in more detail in Chapter 4. 

From the off-axis and laminate test comparisons, it is obvious that, overall, parameter set 2 
correlates best with experimental results. With this parameter set chosen, a complete set of 
endochronic constants has been derived. The complete set of endochronic parameters are 
shown in Table 3. From this point on, these are the parameters used in all analyses. 


2.5 Summary 


This chapter has shown the development of a full set of three-dimensional endochronic 
equations, the development of the incremental form of these equations, and the experimental 
determination of a complete set of endochronic constants. The endochronic equations, based 
on the Conservation of Energy and the Dissipation Inequality, are slightly different in form 
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Table 3. "Best Fit" Set of Endochronic Parameters for AS4/3502 Graphite-Epoxy 



— 

5.12687E-08 

D 

= 

-1.70722E-08 


— 

-4.3511 IE-07 

B 

= 

6.55133E-07 

D 

=s 

2.17931E-06 


= 

1.18373E-06 

Sn 

= 

1.17136E-14 


= 

1.22990E-13 

s 23 

= 

3.28159E-12 

S 2 2 

= 

3.28159E-12 

D 

= 

O.OOOOOE-OO 

B 

= 

3.28159E-12 

Bn 

= 

1.90645E-16 

B 

= 

3.14838E-04 

DO 

D 66 

= 

1.00000 

ni 

= 

2.44791 E-00 

n 2 

= 

2.45967E-00 

n« 

= 

3.18540E-00 


where 


Ajj units are 1/psi 

All other parameters are dimensionless 
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from Pindera and Herakovich's 19 equations, due to a softening rather than stiffening behavior 
for compression loading in the fiber direction. 

The incremental form of the endochronic equations allows them to be used in a finite element 
program. During each increment of loading, stress components are assumed to be propor- 
tional to one another and to the loading increment. This allows explicit integration of the 
equations. The final incremental form of the equations is a summation expression which may 
be implemented into a finite element or laminate analysis program. 

Determination of the endochronic constants required analyzing uniaxial compression tests, 
off-axis tension tests, and losipescu shear tests. The final set of eighteen endochronic pa- 
rameters obtained is unique to AS4/3502 graphite-epoxy. These constants are further unique 
in that they were derived assuming compression loading in the fiber direction. 
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3.0 Finite Element Program AMNISAP 


3.1 Introduction 


Out-of-plane interlaminar stresses are present in boundary layer regions near all free edges 
of notched angle-ply laminates subjected to in-plane loading. In order to account for this in 
a finite element analysis, three-dimensional elements are required. This study uses a three- 
dimensional, material nonlinear finite element program which is a modification of an existing 
linear elastic three-dimensional program \ The present program is called AMNISAP, which 
is an acronym for Anisotropic Material Nonlinear Inversion Symmetric Analysis Program. 

Because of the extreme stress gradients present in the immediate region of the hole, a fine 
mesh was used in this region. In order to handle the large storage requirements produced 
by the fine mesh, a large capacity computer was necessary. Storage and run time efficiency 
were increased by utilizing an inversion symmetry condition and by running the program on 
a computer with vector processing capabilities. 
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Incorporation of the endochronic constitutive equations into a finite element program requires 
solving a set of nonlinear equations. There are many algorithms which may be used to solve 
such equations; each has its own advantages and disadvantages. This chapter examines 
some of the commonly used algorithms. The two algorithms used in the finite element pro- 
gram AMNISAP are developed and described in detail in this chapter. 


3.2 Inversion Symmetry 


Inversion symmetry is a symmetry condition present in all anisotropic shells and plates whose 
geometry, elastic coefficients, boundary conditions and loading are 180 degrees rotationally 
symmetric about some axis 45 . This axis of symmetry, commonly referred to as the center of 
symmetry, is normal to the body and intersects the origin of the in-plane coordinate axes. The 
two-dimensional inversion symmetry transformation consists of rotating the shell or plate 
through 180 degrees about the axis of symmetry. Defining an x, coordinate system (i = 1 ,2,3), 
with the x 3 coordinate the axis of symmetry, a symmetry transformation changes the coordi- 
nates x a (a = 1,2) of each material point of the shell or plate into — x 0 . Figure 22 shows an 
anisotropic shell exhibiting the required conditions for inversion symmetry. 

A plate's requirements are less stringent than a shell's requirements and are listed below. 

INVERSION SYMMETRY PLATE REQUIREMENTS: 

GEOMETRY 


d(x 1 ,x 2 ) = d( -x 1t -x 2 ) 


ELASTIC COEFFICIENTS 
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^qPyp^I’ x 2^ ^aPyp( x 1> x 2^ 


^apvp( x 1> x 2^ P apYp( x 1' x 2^ 


^aPYp^ x 1 ,x 2) ^aPYP^ X 1 > x 2) 


BOUNDARY CONDITIONS 


B(x 1 ,x 2 ) = B( — x v — x 2 ) 


LOADING 


P(x 1f x 2 ) = P( -X 1t -x 2 ) 


M(x 1 ,x 2 ) = M( -x v -x 2 ) 


In the above, d is the distance from the center of symmetry to the plate boundary, C ofJyp are the 
extensional stiffness coefficients, F oPyp are the bending-extension interaction coefficients, and 
D aPyp are the bending stiffness coefficients. B represents the boundary conditions, P repres- 
ents the in-plane loads, and M represents moment (bending) loads. The response of an in- 
version symmetric plate to this loading will follow the relationships below. 


INVERSION SYMMETRIC PLATE RESPONSE: 

DISPLACEMENTS 


U 1 (x.,,x 2 ) = U 1 ( -x 1f -x 2 ) 


U 2 (x-i- x 2 ) = U 2 ( x 1 , -x 2 ) 


W (x 1( x 2 ) = W ( -x.,, — x 2 ) 


STRESS RESULTANTS 
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^ap (*l> *2) ®ap ( *ii *2) 

^o3 (^1 > ^ 2 ) ^a3 ( X 2 ) 

ct 33(*1.X2) = CT 33 ( -X^ ~X 2 ) 
where a, p = 1, 2 




In the above, U, and U 2 represent the in-plane displacements in the Xi and x 2 directions, and 
W represents the out-of-plane displacements. The stresses a op are the in-plane components, 
ct o3 are the out-of-plane shear stress components, and a 33 is the out-of-plane normal stress. 

From these requirements, it is apparent that an angle-ply laminated plate with a central hole 
under uniaxial loading (Figure 23) exhibits inversion symmetry. The center of symmetry co- 
incides with the center of the hole. 

By combining the inversion symmetry condition with the mid-plane symmetry condition which 
is present in all symmetric laminates, only one quarter of the plate need be analyzed 
(Figure 23). The mid-plane symmetry condition is invoked simply by fixing the mid-plane 
nodai points' out-of-plane displacements W. Exploitation of the inversion symmetry condition 
requires modification of the finite element program. Appendix D gives a detailed description 
of the modifications required. It is desireable to add a flag variable for each nodal point of a 
finite element grid identifying those nodes which are independent and those which are de- 
pendent. Once these modifications are made, both inversion symmetric and non-inversion 
symmetric analysis may be performed with the same program. 
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3.3 Vectorization 


Even with the inversion symmetry condition fully utilized, the storage requirements for many 
of the analyses performed in this study (dependent upon mesh refinement) were excessive for 
the IBM 3090 or 3084 main frame computers at Virginia Tech. Therefore, the Cyber 205 super 
computer at the NASA Langley Research Center was used. This machine easily handled the 
storage requirements. An added bonus of the Cyber 205 is its parallel and vector processing 
capabilities. Using vector instructions, one instruction performs many operations, thus re- 
ducing CPU time significantly. The measure used to define this performance feature is 
MEGAFLOPS (Millions of Floating Point Operations Per Second). The Cyber 205 operates in 
the range of 50 to 200 MEGAFLOPS, depending upon the vector length. Appendix E shows an 
example subroutine comparing scalar operations to vector operations. 

Because many do-loops in a program are not candidates for vectorization and because there 
are usually only a few do-loops which use the majority of CPU time, it is necessary to be se- 
lective when vectorizing a program. To determine which scalar code subroutines were the 
most time consuming, a timing option was run which output the total time and time per call 
for each subroutine. Table 4 shows these results. It is noted that these results vary de- 
pending upon job size. Table 4 shows that subroutines COLSOL (the global equation solver) 
and ELKAY (calculates the element stiffness matrices) used 85% of the total CPU time for the 
example considered. It is in these subroutines that the largest gains were realized through 
vectorization. 

Subroutines COLSOL, ELKAY, and STRESS were vectorized by an automatic vectorizing utility 
named VAST. The FORTRAN program was run through this utility before it was compiled. 
VAST replaces scalar do loops with vector syntax wherever possible. It also returns a listing 
file which shows those do loops that were vectorized, lists those do-loops that weren't 
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Table 4. Timing Results for Scalar Linear Code Using a 20 Element Mesh 


Subroutine 


ADDBAN 


ADDRES 


AMNISAP 


ASSEM 


BCIMP 


BEDGE 


BSTRS 


CLEAR 


COLHT 


COLSOL 


DMAT 


DN16 


DN20 


DN24 


ELKAY 


ELSUM 


ERROR 


GCOORD 


GPN16 


GPN20 


GPN24 


INPUT 


INVERT 


LOADS 


LOADV 


STRESS 


SYMDUMP 


WRITE 


Call Freq. 
Counts 

% Call 
Freq. 

Total CPU 
Time (secs) 

CPU Time 
Per Call 

20 

6 

0.052793 

0.002563 

1 

0 

0.000767 

0.000488 

0 

0 

0.015310 

0.000000 

1 

0 

0.045322 

0.031250 

1 

0 

0.017408 

0.015625 

20 

6 

0.000020 

0.95E-06 

1 

0 

0.040216 

0.031250 

1 

0 

0.009775 

0.007813 

20 

6 

0.002568 

0.000122 

2 

0 

1.73258 

0.750000 

2 

0 

0.001442 

0.000488 

2 

0 

0.000342 

0.000122 

2 

0 

0.000002 

0.95E-06 

2 

0 

0.000002 

0.95e-06 

20 

6 

3.78203 

0.187500 

2 

0 

0.042023 

0.015625 

3 

1 

0.000006 

0.000002 

4 

1 

0.000060 

0.000013 

2 

0 

0.000159 

0.000061 

2 

0 

0.000002 

0.95e-06 

2 

0 

0.000004 

0.000002 

i 

0 

0.174491 

0.125000 

180 

53 

0.004861 

0.000027 

1 

0 

0.000744 

0.000488 

1 

0 

0.002422 

0.001953 

1 

0 

0.48726 

0.250000 

1 

0 

0.043696 

0.031250 

44 

13 

0.001848 

0.000042 

1 

0 

0.046071 

0.031250 


% Total 
CPU Time 


TOTAL CALLS = 340 


TOTAL CPU TIME = 6.499723 secs 
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vectorized, and explains why certain do-loops were not vectorized. In many instances, simply 
rearranging a do-loop will allow VAST to vectorize it. In other cases, the algorithm must be 
rethought and changed. Besides using VAST on the above subroutines, a few do-loops 
throughout the program were vectorized explicitly using vector syntax. The major gains, 
however, were obtained through using VAST. 


3.4 Program Efficiency Comparison 


Table 5 shows vectorized timing results for a 20 element mesh (very coarse) model of a 
compression loaded angle-ply plate with a central hole. The comparison is for linear elastic 
analysis with no inversion symmetry. It must be noted that significantly larger speed factors 
are possible using a more refined mesh due to larger array sizes. Table 5 shows an overall 
speed factor of almost 4 with very little effort. Most of these gains were realized by rear- 
ranging do-loops within the subroutines. 

Table 6 compares inversion symmetry analysis (which requires a quarter plate) with half plate 
and fixed node quarter plate analysis for a compression loaded [ ± 30], laminate with a 
central hole. The fixed node quarter plate analysis is performed with the nodes along the y- 
axis fixed in the x-direction. These are approximate boundary conditions which give errone- 
ous results in regions close to the y-axis but give acceptable results elsewhere \ Inversion 
symmetry results agree with half plate results everywhere. Note that the half plate storage 
requirements are more than 2.5 times the inversion symmetry requirements, while the inver- 
sion symmetry requirements are only 1.2 times greater than the fixed node quarter plate re- 
quirements. Half plate CPU time is more than 3.0 times that of inversion symmetry, yet 
inversion symmetry CPU time is only 1.3 times greater than that of the fixed node quarter plate 
analysis. Half plate total run cost is 2.9 times that of inversion symmetry, while inversion 
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Table 5. Vectorizing Timing Results for Non-Inversion Symmetric Linear Elastic Analysis Using a 
20 Element Mesh 


Vectorization 

Modifications 

Total 
Run Time 
(CPU'S) 

Overall 

Speed 

Factor 

COLSOL 

Speed 

Factor 

ELKAY 

Speed 

Factor 

STRESS 

Speed 

Factor 

AMNISAP in original form - 
no vectorization 

6.50 

1.00 

1.00 

1.00 

1.00 

VAST applied to AMNISAP 
(only subroutines ELKAY, 
COLSOL, STRESS) 

3.49 

1.86 

5.44 

1.65 

1.22 

VAST applied to AMNISAP - 
2 do-loops in ELKAY 
modified (80 & 85) 

3.18 

2.04 

5.42 

1.91 

1.22 

VAST applied to AMNISAP - 
2 more do-loops in ELKAY 
modified (60 & 70) 

3.02 

2.15 

5.39 

2.09 

1.22 

VAST applied to AMNISAP - 
2 more do-loops in ELKAY 
modified (90 & 91) 

1.75 

3.72 

5.41 

6.97 

1.22 

VAST applied to AMNISAP - 
do-loop 100 modified 
in COLSOL 

1.66 

3.93 

7.61 

6.97 

1.22 


Speed Factor = . ° g!l R™ Time 
Vectorized Run Time 


Finite Element Program AMNiSAP 


80 





































Table 6. Inversion Symmetry Results for [±30], Laminate with a Central Hole Run on Scalar 
Code 


Parameters 

Half 

Plate 

Inversion 

Symmetry 

Fixed Node 
1/4 Plate 

# of Elements 

200 

100 

100 

# of Nodes 

2004 

1041 

1041 

# of Equations 

5338 

2669 

2698 

# of Global Stiffness Terms 

1,809,178 

630,897 

406,513 

Total Storage (words)* 

2,886,800 

1,163,400 

943,400 

Maximum Half Bandwidth 

5320 

2606 

366 

Mean Half Bandwidth 

338 

236 

150 

CPU Time* 

324.1 

104.8 

71.4 

Run Costs (CRUs)* 

356.4 

124.3 

99.2 


* NASA Langley's Cyber 205 
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symmetry analysis' total run cost is only 1.3 times that of fixed node quarter plate analysis. 
In conclusion, for a fraction more of the fixed node quarter plate computer resources, half 
plate results are realized by using inversion symmetry at a significant reduction in computer 
resources. 

If the results from Table 5 and Table 6 are combined, a total speed factor of 12 is realized 
when comparing half plate non-vectorized analysis to inversion symmetric vectorized analy- 
sis. Note also that this result is conservative being based on the speed factor for the 20 ele- 
ment mesh. 


3.5 Numerical Solution of Nonlinear Algebraic Equations 


Nonlinearities in structural mechanics exist in two forms - material nonlinearities and ge- 
ometric nonlinearities. Material nonlinearities arise from nonlinear constitutive laws. Ge- 
ometric nonlinearities are the result of finite changes in the geometry of the body. Many 
numerical solution techniques are currently available for solving the set of nonlinear algebraic 
equations arising from both causes. However, all of these techniques have limitations. Some 
methods work well for geometric nonlinearities while others work best for material nonline- 
arities. The accuracy and efficiency of each method varies widely as well. All of these factors 
and features must be accounted for when choosing a method for a particular application. 

This study focuses on material nonlinear finite element analysis. The solution techniques 
most commonly used for material nonlinearities are: (1) Incremental or stepwise procedures, 
(2) Iterative procedures, and (3) Step-iterative or mixed procedures. Each of these methods 
has its own distinct advantages and disadvantages and each one was considered as a sol- 
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ution technique in this study. Common to all of these techniques is the set of equations they 
must solve. This set is typically written for a single element as 

[K]{U} = {F} (3.1) 


where 

[K ] = [K] + [K NL ] 

[k nl ] = [K NL ({U> , {F»] 

{11} = nodal displacements 
{F} s nodal forces 
Ck] a stiffness matrix 
[K] a linear portion of stiffness matrix 
Ck nl : a nonlinear portion of stiffness matrix 

Figure 24 shows a typical force/displacement curve for a given element and a typical nonlin- 
ear stress/strain curve. Following is a discussion of the solution techniques. 


3.5.1 Incremental Procedures 


Incremental procedures involve dividing the load into a series of increments. A solution is 
generated for each load increment until the final load is reached. Normally, the load incre- 
ments are equal in magnitude but this is not a requirement. The stiffness matrix [R] is a 
constant for a given load step, but may change for each load increment. This matrix is re- 
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(a) Force-Displacement Curve 



(b) Stress-Strain Curve 


Figure 24. Typical Nonlinear Curves 
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ferred to as the tangent stiffness matrix in this procedure. For each incremental set of loads, 
a set of incremental displacements are generated. The final set of loads and displacements 
is the sum of the incremental loads and displacements. This is written symbolically as 


{F} = {F 0 } + X {AF,> 
1 = 1 


(3.2) 


M 


{U> = {U 0 } + I {AU,} 


where the A symbol represents an increment, M is the total number of increments, and sub- 
scripts 0 refer to the initial state. At a given load step, 


{Fj> = {F 0 } + j S {AF,} 


(3.3) 


{Uj> - {U 0 > + S {AU,} 
! = 1 


The incremental element equations are then written as 

[K,-i] {AU,} = {AF,} for i- 1.2.3 M (3.4) 

where 

[K|-l] = [K i -i({U,_ 1 > . {F,.,})] 

Figure 25 depicts this incremental procedure. 

The accuracy of the incremental method is very dependent upon the load step size. The ac- 
curacy can be improved by taking smaller load steps. However, as load step size decreases 
and the number of load steps increases, the efficiency of the procedure decreases because 
[K,_J must be computed for each load step. 
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3.5.2 Iterative Procedures 


Iterative procedures consist of a sequence of calculations in which force equilibrium is more 
exactly satisfied with each new iteration. The sequence of calculations (or iterations) is con- 
tinued until the force imbalance (or residual) falls below an acceptable magnitude. The full 
force is applied at each iteration of this method. After each iteration, the force residual is 
calculated and used in the next iteration to compute the displacement. 

The two most common iterative techniques used in finite element analysis are 46 : (1) Func- 
tional iterations (or successive substitutions) and (2) the Newton Raphson method. In both of 
these methods, the finite element equations are written 

[K]{U} = {F} - {F NL > (3.5) 


where 


[K] = linear portion of stiffness matrix 
{F} = applied load vector 

{f nl } = psuedo load vector resulting from [k nl ] {U} 

3.5. 2.1 Method of Functional Iterations 

This method seeks to find {U} for a given {F }. It starts with an initial estimate of the dis- 
placements {U> (solving the linear part of the finite element equations yields a reasonable 
estimate) and performs the following iteration sequence 

[K]{u< r >} = { F> - (F NL ({u (r - 1 >>)} (3 .6) 
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where r, r - 1 refer to the iteration cycle. The sequence of iterations is continued until the 
difference in displacements between successive iterations falls below a preset value. The 
linear stiffness matrix [K ] is decomposed only once and remains constant throughout the it- 
erative process -- the changing nonlinear terms are contained in the psuedo load vector 
{F NL >. Figure 26 shows a schematic representation of this iterative procedure for a one degree 
of freedom system. 

3.5 .2.2 Newton Raphson Method 

In the Newton Raphson method, a force residual represents the initial unbalance between the 
two sides of equation (3.5). This force residual is written 

{f(U)> = [K]{U} - {F} + {F NL > (3.7) 

This method expands the force residual (f(U) } in a Taylor Series about the initial estimate. 

{f(U + AU)> = <f(U)> + + Error ({AU}) 2 (3.8) 

An approximation used in the Newton Raphson Method is to truncate the series after two 
terms which yields the error term in the above expansion. A requirement of the method is to 
set the expanded force residual to zero, {f(U + AU)} = {0} . These approximations yield 

t^TT-KAU} - -{f(U)> (3.9) 

d{U> 

Now from equation (3.7), with {F} independent of {U} 

= IK] + = IK] + IK(U)] (3.10) 
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Combining (3.9) and (3.10) gives 


[K + K(U)] {AU> = -{f(U)> 

In terms of the iterative cycle, this is written 

[ [K] + [K(u (r-1) ) ] ]{AU (r) > = - {f(u (r-1) )} (3.11) 

Note that each iterative cycle of the Newton Raphson Method requires computation of {f(U)} 
, computation of [ K(U) ], and decomposition, forward reduction, and back substitution to 
compute { AU }. Figure 27 shows a schematic representation of this iteration procedure, 
where {K r } is the tangent stiffness matrix and the superscript r represents the iteration num- 
ber. 


3.5.3 Step-Iterative Procedures 


Step-iterative procedures are simply a combination of incremental procedures and iterative 
procedures. These method yields higher accuracy but may require more computational effort. 
A typical step-iterative method is depicted in Figure 28. More details on this method are 
presented in a later section. 


3.5.4 Procedure Comparison 


The advantages of the incremental procedures are many fold, not the least being simplicity. 
They are easy to program, they admit a simple physical interpretation, they are applicable to 
nearly all types of nonlinear behavior, and they provide a relatively complete description of 
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Figure 28. Force/Displacement Curve for a Typical Step-Iterative Procedure 


Finite Element Program AMNISAP 


92 


the load-deformation behavior. However, their disadvantages are just as numerous. Incre- 
mental procedures require solving a linear system of algebraic equations at each step which 
necessitates formation of [K], and requires decomposition, forward reduction, and back sub- 
stitution to compute {AU >. Also, there is no rational way of selecting step size or checking 
convergence of the solution; hence the accuracy of the solution tends to drift. 

The advantages and disadvantages of the iterative methods must be separated to include 
those of the functional iterative method and those of the Newton Raphson method. The ad- 
vantages of the functional iterative method are its simplicity of application -- the stiffness 
matrix [K ] is decomposed only once - and the solution's accuracy may be specified through 
the convergence criterion. However, the functional iterative method is not suitable for highly 
nonlinear problems. It displays slow convergence (and sometimes divergence) for elasto- 
plastic or near elasto-plastic problems, and the solution converges in a linear fashion requir- 
ing many iterations. 

The Newton Raphson method's advantages are that it works well for highly nonlinear prob- 
lems, and that it converges quadratically requiring fewer iterations (not necessarily faster CPU 
time). Its disadvantages include slow convergence for highly nonlinear problems, conver- 
gence problems due to oscillation about a local maximum or minimum point, and a high 
computational expense because [K + K ] is decomposed many times. 

The step iterative procedures combine the advantages of the incremental and iterative meth- 
ods and tend to minimize the disadvantages. Because of this, step iterative methods have 
become increasingly prevalent in modern nonlinear finite element analysis. 


Many factors were considered when choosing a nonlinear solution technique for the finite el- 
ement program of this study. The first of these factors was the complexity of the endochronic 
equations. Explicit differentiation with respect to the displacements, as required in the Newton 
Raphson method, is very difficult if not impossible with the form of the endochronic equations. 
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Numerical differentiation requires still more complexity and computational time as well as an 
additional approximation in the solution. The form of the endochronic equations was also 
considered - strain is written in terms of stress. Although a rederivation of the equations 
could yield stress in terms of strain, the endochronic parameters are much more difficult to 
determine with the equations in this form. Also, as shown in Chapter 2, the equations can only 
be integrated incrementally, with proportional loading assumed over each increment. There- 
fore, the form of the endochronic equations limited the choice of solution techniques. The final 
factor considered was the size of the finite element analysis. A method requiring many de- 
compositions, forward reductions, and back substitution of the stiffness matrix would consid- 
erably increase CPU time. 

With all these factors considered, two step-iterative procedures were chosen — the constant 
stress method and the constant strain method. These methods are discussed and developed 
below. 


3.5.5 Constant Stress Method 


The constant stress method is an initial strain method in that the total strain is separated into 
a conservative portion and a non-conservative portion. The non-conservative strain (also 
called the irreversible strain) is integrated separately as a psuedo force array which is added 
to the applied force array to yield the total force array. Therefore, all the nonlinearities in the 
equations are represented by this psuedo force array. Because of this, the stiffness matrix 
IK ] is decomposed only once in the solution process. This reduces CPU time and increases 
computational efficiency. The initial strain formulation of the finite element equations is shown 
below. 


{e > = {e R } + {e 'R> 


(3.12) 
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where 


{s R > = [A]{a> 
{e IR } = {e IR (z)> 


{e } - reversible portion of strain 
{e IR > = irreversible portion of strain 
z = intrinsic time scale 

Therefore, 

{a} = (A] 1 {e R > = [A] _1 {e - e IR } (3.13) 

The strain energy of an element (A*) is expressed as 

A6 = Jv4 {eR}T{a > dV = Jvy({">V> - {e ,R > T {0>)dV (3.14) 

where 

v e = volume of an element 

Combining (3.13) and (3.14) gives the strain energy in the form 

A ' - Jv^jWVr’M - -i{e} T [A]- , {e"'> - J-WlAJ-Vt) + •l{e IR }[A]- V R }) dV 
= jv-y(< e > T l A r'(E> - 2{e} T |AJ -1 {e ,R > + {t'TlAr ’{>:"*}) dV (3.15) 
The strain-displacement equations can be written 
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{e} = [B]{U> 


{3.16) 


where [B] is a matrix containing the derivatives of the finite element interpolation functions. 
The transpose of the strain vector is 

{<0 T = {U} T [B] t 


Therefore, 


A e - J v »-y({ u } T tB] T [Al ^BKU) - 2{U} T [B] T [A]" 1 {e ,R > + {e IR > T [A]' 1 {e IR }) dV (3.17) 


The potential energy of an element due to externally applied forces is given by 

T e = -{U} T {P e > 


(3.18) 


where 


{P e } = applied load vector for a given element 

The total potential energy of an element (11*) is the sum of the strain energy (A*) and potential 
energy (r*). 


n e = A® + r e (3.19) 

The total potential energy of the finite element model is the sum of the elemental potential 
energy over the elements. 

N 

n = in® (3.20) 

e = 1 


Therefore, the total potential energy is 
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" " j,[ J'vl< U > T Pl T l A l' , |BKU}dV - j^Ww-'Wd* ] + (321| 

Jj J v .{e l "»r , {t ,R ><IV - <U> T <P'> 1 


where 


N = number of elements 


In order to obtain the correct displacement solution solution from (3.21), the total potential 
energy IT must be minimized. This is done by taking the partial derivative of the total potential 
energy with respect to the displacements as follows 


^fj}" " {0} " J^PlWlBl dV {U> - dV - {P e } ] 


(3.22) 


Defining: 


[K e ] = J v .[B) t [A]" 1 [B] dV 


and 


{l" 1 } = j^BlV] 1 {e IR } dV 

{(•> . {f™-} + <p« } 


where 


[K e ] = element stiffness matrix 




element psuedo load vector 
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{p e } = element force vector 


{t®} = element total load vector 


gives 


I {[K e ]{U} - {f}} = {0} 
e = 1 


where 


N = number of elements 


Writing the above in global matrix form yields 


IK]{U> = {F> 


where 


{F> = {F NL > + {P> 


N 

{f nl } = £ f ,[B] t [A] _1 {e IR } dV s global psuedo force array 

e = 1 JV 


IN 

{P} = L {P e } = global applied force array 


e= 1 


The iterative form of equation (3.23) is written as 

IK] {U< r) } = { FNL ({Auf r_1) })} + {P|> 


where 


r = iteration cycle 


(3.23) 


(3.24) 
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i = load step number 

Note that the full load is applied at each step resulting in total (as opposed to incremental) 
displacements. This differs somewhat from true incremental procedures in which only an in- 
cremental load is applied at each step. Although the full load and full displacements are re- 
presented in the equations, the psuedo force array {F NL } is a function of the incremental 
stresses {Act } at the previous iteration. 

With the equations written in the above form, the constant stress method simply defines the 
order of calculations. This order of calculations is outlined below: 

1) For a given load step i and iteration r, convergence is checked comparing the total dis- 
placements at the present iteration {U| rl } with the total displacements at the previous iteration 

{u M>> . 

2) The element strains are computed from the element displacements 

{E, (r) > = (B]{u[ r) } (3.25) 

3) The stresses in the element are found from the strains at the current iteration {e! r) } and the 
irreversible strains at the previous iteration {eS R,r-1) } 

{o| r| }-|A1-'{{ E r'}-{e,' R ''-"}} (3.26) 

4) The incremental stresses in the element are calculated by subtracting the stresses from the 
previous load step from the stresses at the present load step i 

{Aa| r! } = {CT, (r) > - {CTp^} (3.27) 
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5) The irreversible strains from the present iteration r are computed using the incremental 
stresses {Aof'} 

{e! R|r| }- {^({Ac!' 1 })} (3.28) 

6) Finally, the psuedo forces for the r + 1 iteration are calculated from the irreversible strains 
at the present iteration r 

{ F| NL(r + 1)j = {F NL ({E! R(r) })J. (3.29) 

If displacement convergence is satisfied, the load is incremented and the process is repeated. 
If displacement convergence is not satisfied, the above sequence is repeated in a do-loop until 
a convergence criterion on the displacements is satisfied. The constant stress method gets 
its name because the sequence of calculations are performed at essentially a constant stress 
level. The displacements and thus the strains change during the iteration process until con- 
vergence is satisfied. Figure 29 gives a schematic representation of the constant stress 
method and Figure 30 shows a flow chart for the method. 


3.5.6 The Constant Strain Method 


The constant strain method is also an initial strain method and uses the same form of the finite 
element equations as the constant stress method. It differs only in the sequence of calcu- 
lations required to solve the finite element equations. Whereas the constant stress method 
iterations were performed at essentially a constant stress level, the constant strain method 
iterations are performed at essentially a constant strain level. The stresses and thus the ap- 
plied loads change during the iteration process until a convergence criterion is satisfied. 
Whereas the constant stress method consists of one do-loop, the constant strain method 
contains two do-loops - an inner loop and an outer loop. The inner loop iterates at a given 
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Calculate and assemble element stiffness 
matrices and nodal forces 
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Figure 30. Flow Chart Showing Constant Stress Method 
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strain level until a stress convergence criterion is satisfied. The outer loop iterates until a 
displacement convergence criterion is satisfied. Once the the outer loop has converged, the 
procedure is incremented. The constant strain method is outlined below: 

1) For a given load step i and iteration r, convergence is checked comparing the total dis- 
placements at the present iteration {U| rl } with the total displacements at the previous iteration 

{ ur’> }. 


2) The element strains are computed from the element displacements 

{ 4 °} - |BJ{U«} 


(3.30) 


3) The stresses in the element are calculated from the strains of the current iteration {ej r) } and 
the irreversible strains at the previous iteration {el"*'- 1 '}. 

{a!' 1 } = W’{K’} - {=!""*"}} (3.31) 

4) Convergence is checked comparing the total stresses at the present iteration {a| rl } with the 
stresses at the previous iteration {o| r-1) }. 

5) The incremental stresses in the element are calculated by subtracting the stresses from the 
previous load step from the stresses at the present load step i 

{Aof r) } = {o[ r) > - {oj'ij (3.32) 

6) The irreversible strains for the present iteration r are computed using the incremental 
stresses {Aaj r) > 

{e| R(r) }= ^ E IR ({Aa{ r) })^ (3.33) 
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7) If stress convergence is satisfied, then the solution proceeds. If it is not satisfied, the sol- 
ution returns to step 3. 

8) The psuedo force array for the r + 1 iteration is calculated from the irreversible strains at 
the present iteration r 

^pNLfr + IJj. ^F NL ({Ej' R(r, })J. (3.34) 

If displacement convergence is satisfied, the load is incremented and the process is repeated. 
If displacement convergence is not satisfied, the equations are solved for the r + 1 iteration 
and the process is repeated. Figure 31 illustrates a schematic representation of the constant 
strain method and Figure 32 shows a flow chart. 


3.6 Summary 


The three-dimensional material nonlinear finite element program AMNISAP contains both the 
constant stress and the constant strain methods of solution. The user may choose the desired 
method. The reason for incorporating both methods into the program is that the constant 
stress method works better for some problems while the constant strain method works better 
for others. 

The constant stress method has been previously shown to be inherently unstable under cer- 
tain circumstances 46 ; it may suddenly and catastrophically diverge. However, under stable 
conditions, it converges rapidly. On the other hand, the constant strain method is inherently 
stable, but converges much slower than the constant stress method. It too may diverge, but 
does so in a progressive, not catastrophic, manner. 
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Figure 31. 


Force/Displacement Curve for Constant Strain Method 
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Figure 32. Flow Chart Showing Constant Strain Method 


Finite Element Program AMNISAP 


Iteration Limit Exceeded 





AMNISAP is a multipurpose program with options to run inversion symmetric or non- inversion 
symmetric analysis, and linear material or nonlinear material analysis. Within the nonlinear 
material option, the user may choose between two solution methods -- the constant stress 
method or the constant strain method. Finally, the program is written in Fortran 77 so that it 
is interchangeable with any computer facility having a Fortran 77 compiler. 
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4.0 Laminate Response 


4.1 Introduction 


This chapter presents an analysis of unnotched and notched angle-ply laminates loaded under 
uniaxial compression. In each analysis, analytical and experimental results are compared. 
Also, where appropriate, comparison is made between material linear and material non-linear 
predictions. The angle-ply laminates considered experimentally are: [( ± 10) 12 ], , 

[( ± 20) 12 ], , [( ± 30) 12 ] s , and [( ± 45) 12 ], . Although these test laminates had 48 plies, they 
were modeled in the finite element analysis as having only 4 plies {i.e. [ ± 10], , [ ± 20], , 
[ ± 30], , and [ ± 45], ) due to computer storage limitations. The analysis of unnotched 
laminates does not include edge effects and gives identical classical lamination theory results 
whether 48 plies or 4 plies are modeled. On the other hand, the analysis of laminates with 
holes does include edge effects around the hole boundary. These edge effects are dependent 
on the laminate stacking sequence. Since 48 ply laminates could not be investigated, it was 
decided to investigate four ply laminates which have the same alternating ± <p sequence as 
the 48 ply laminates. It is believed that the overall trends in the 4 ply laminates are similar 
to those in the 48 ply laminates. It is also important to note that the hole diameter is not the 
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same for ail the notched laminates. Table 7 gives the dimensions of the laminates tested. 
The dimensions of the finite element models are identical to those in Table 7 except for the 
thickness dimension which was 0.01" in all finite element models. 

Note that the [( ± 10) u ], and [( ± 20) 12 ], laminates with 1/2" diameter holes have a w/d 
(width to diameter) ratio of 5 while the other laminates have a w/d ratio of 10. The w/d ratio 
of 5 was used to investigate the influence of the width to diameter ratio on failure. In order 
to differentiate between the unnotched and notched laminates, they are referred to as 
[ ± <p], and [ ± (<p)i 2 ]„ respectively. 


4.2 Unloading Considerations 


In the course of performing the finite element analysis, it was observed that the solution would 
converge quite rapidly for each load step up to a critical load level; however, at this critical 
load level, the solution would diverge rapidly. For some cases, this divergence phenomenon 
would occur before the load reached the experimentally observed failure load. Close inves- 
tigation into this problem revealed that a component (or possibly many components) of stress 
was unloading at the critical load. Because the unloading causes a change in the sign of the 
slope (Figure 33) on the component's stress-strain diagram, the solution algorithm (either 
constant stress or constant strain) can not find the correct value of stress. The solution os- 
cillates back and forth (with each oscillation growing larger) in an attempt to find the correct 
Act,. Figure 33 shows this phenonomen schematically for the constant strain algorithm. 

In Figure 33, the symbols oC, and c, represent the stress level at the critical load and after the 
critical load step, respectively. The letter i represents the iteration number. At the critical 
load step, eC represents the strain from the initially applied load, and because the irreversible 
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Table 7. Test Specimen Dimensions 


Layup 


# of Hole 

Specimens Length Width Thickness Diameter 

Tested (in.) (in.) (in.) (in.) 



l ± 10]. 

2 

10.0 

2.0 

0.25 

t ± 20], 

2 

10.0 

2.0 

0.25 

I ± 30], 

1 

10.0 

2.0 

0.25 

[ ± 45], 

1 

10.0 

2.0 

0.25 




l ( ± 10)„], 

2 

10.0 

2.5 

0.25 

0.25 

K ± 10)*]. 

2 

10.0 

2.5 

0.25 

0.50 

l ( ± 20) 12 ], 

2 

10.0 

2.5 

0.25 

0.25 

[( ± 20), 2 ], 

2 

10.0 

' 

2.5 

0.25 

0.50 

[( ± 30) 12 ], 

1 

10.0 

5.0 

0.25 

0.50 

[( ± 45),*]. 

1 

10.0 

5.0 

0.25 

0.50 
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strain has not yet been calculated, ctC, corresponds to the linear elastic stress. The irre- 
versible strain eC| R is then calculated from ctC, (3.29). Psuedo loads are calculated from eC( R 
(3.30) and the process is repeated until a convergence criterion is satisfied at the fourth iter- 
ation. The stress oC 4 corresponds to the strain tC. At the following load step, the strain e is 
calculated from the applied loads and psuedo loads. The stress a, is calculated from the 
strain e and eQ R from the previous load step. The irreversible strain is calculated from 
ct 3 , but because the stress-strain curve slopes downward below stress level <x,, the program 
calculates a large indeterminate value for e!, r (imagine the curve continuing with no inflection 
point). At the second iteration, e^ r is used to calculate ct 2 . The irreversible strain e^ r calculated 
from o 2 is very small. At the third iteration, ct 3 is larger than ct, due to the small value of z'f. 
The irreversible strain e^ r calculated from o 3 is larger still than ef and causes the stress ct 4 at 
the fourth iteration to be negative. When the stress becomes negative, divergence proceeds 
rapidly as the upwards and downwards arrows depict. 

Although it is unknown whether this unloading phenomonen actually occurs in these lami- 
nates, the endochronic theory does predict it. Zinov'ev and Sarbaev 21 show a similar 
occurance for an unspecified fiberglass [ ± 40], laminate loaded in uniaxial tension. Figure 
5 in their paper shows the normal stress transverse to the fibers (ct 2 ) to be positive initially, 
but become negative approximately midway through the loading history. Because Zinov'ev 
and Sarbaev used a different formulation of the endochronic equations (i.e. they used multiple 
time scales whereas this study only uses a single time scale), it is highly unlikely that this 
phenomonen is an error in the formulation of the equations or in the finite element program. 
Rather, it appears that it is inherent to the endochronic theory. To prove or disprove whether 
the unloading is a real occurrence would require experimental verification. 

Because the solution divergence phenomenon sometimes occurs before experimentally ob- 
served failure, a condition was imposed in the finite element program which allows the sol- 
ution to progress past the critical load level. For the stress components which unload, the 
stress increments are set to zero (Act, = 0) at the present load step and at all future load 
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steps. This essentially forces the components' stress-strain moduli to zero for the remainder 
of the loading history. In this manner, the strain components are calculated, but the stress 
components are artificially regulated. Of course, imposing this condition is equivalent to 
specifying an unknown, which hinders the solution's accuracy. However, this was found to 
be the most effective way to reach higher loading levels. In future sections, this condition is 
referred to as the zero-slope condition. 


4.3 Experimental Results 


Burns et al 1 performed experiments on the laminates listed in Table 7. The specimens were 
compression loaded quasi-statically in a rigid frame which provided clamped conditions at the 
loaded ends and simply supported conditions at the sides. Two specimens of each laminate 
configuration were tested except where noted in Table 7. Further detail of the experimental 
procedures is outlined in reference 


4.4 Finite Element Analysis 


The material non-linear, three-dimensional finite element program AMNISAP was used to 
generate all the results. The elements used are 16 node, 3-D isoparametric displacement 
based elements. For further details of AMNISAP, refer to Chapter 3. 

The finite element grid used to model the unnotched angle-ply laminates is an inversion 
symmetric 4 element mesh with one element through the thickness of each ply. This coarse 
mesh is not sensitive to edge effects and gives results identical to classical lamination theory. 
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The finite element grid used to model the notched angle-ply laminates with 1/2" diameter 
holes is an inversion symmetric 744 element mesh; the 1/4" diameter hole laminates were 
modeled with an inversion symmetric 792 element mesh. The meshes for the notched lami- 
nates have three elements through the thickness of each ply giving a total of six elements 
throughout the thickness of the laminates. 

The ends of the laminates were loaded with a force distribution that gave each of the end 
nodes an approximately equal displacement. The end loading was modeled with end forces 
rather than end displacements so that it was not necessary to decompose the stiffness matrix 
for each iteration of the solution. 

The load step size for each of the laminates was based upon the experimental failure loads. 
Twenty load steps were used as a compromise between accuracy and computational effi- 
ciency. The load step size for each laminate was determined by multiplying the experimental 
failure loads by approximately 1.2 and dividing by 20. The factor of 1.2 was used to allow the 
analyses to load past the experimental failure load if the solution proceeded for 20 load steps. 
An exception to this was used in analyzing the [ ± 45], unnotched laminate. The results for 
this laminate exhibited more load step size sensitivity than the other laminates. In order to 
best match the experimental global stress-strain results, only 7 load steps were used. 
Table 8 shows the number of elements, number of nodal points, number of load steps, far field 
stress per load step, and the displacement based convergence criterion magnitude for each 
of the laminates analyzed. 

The constant stress solution algorithm was used in all of the runs. This algorithm converged 
much faster than the constant strain algorithm. For laminates whose solution did not proceed 
to the experimentally observed failure load, the constant strain algorithm was used in an effort 
to advance the solution further. In all cases, no improvement was observed. 
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Table 8. Finite Element Parameters for Angle-Ply Laminates with and without Holes 



Layup 

Hole 

Size 

#of 

Elements 

# of 
Nodes 

#of 

Load Steps 

F.F.S. * 
per Load 
Step (psi) 

Displacement 

Convergence 

Criterion 

• 

{ ± 10J. 

- 

4 

39 

20 

5760 

0.001 


[ ± 20], 

- 

4 

39 

20 

2688 

0.001 

• 

[ ± 30]. 

- 

4 

39 

20 

1920 

0.001 


[ ± 45]. 

- 

4 

39 

7 

6067 

0.001 

• 

[( ± io)„]. 

1/4' 

792 

3171 

10 

4608 

0.001 


K ± 20) 1t ], 

1/4" 

792 

3171 

13 

2688 

0.001 

• 









K ± io ) 12 ]. 

1/2" 

744 

2947 

12 

3840 

0.001 

• 

[( ± 20) 12 ], 

1/2" 

744 

2947 

15 

2304 

0.001 


({ ± 30) 12 ], 

1/2" 

I 

744 

2947 

14 

1920 

0.001 

• 

K ± 45) 12 ]. 

1/2" 

744 

2947 

13 

1920 

0.001 


F.F.S. - Far Field Stress 
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4.5 Unnotched Angle-Ply Laminates 


4.5.1 Global Stress-Strain Response 

4.5.7. 7 Axial Response 

Figure 34 shows comparisons of theory and experiment for the axial response of the four 
laminates considered. The experimental curves in Figure 34 are the same as those shown 
in Figure 18; however, the finite element results are given with smaller load steps in 
Figure 34. The curves are included in one figure for ease of comparison. As indicated in the 
figure, the comparison between theory and experiment is quite good. The response of the 
C ± 10].. [ ± 20],, and [ ± 30], laminates is nearly linear to failure whereas the [ ± 45], 
laminate exhibits significant nonlinearity. 

Comparison between theory and experiment is generally good; however, some interesting 
phenomena can be noted. Near the point of failure, the [ ± 10], laminate exhibits more axial 
softening response than the [ ± 20], laminate. This is due in part to the high compressive 
stress a, in the [ ± 10], laminate which causes some fiber direction softening. In addition, 
the [ ± 10], laminate s response is more dependent upon the fiber direction response. The 
fiber direction softening response can be seen in Figure 2 of Chapter 2. Note that at equal 
stress levels, the [ ± 20], laminate exhibits more axial softening response than the [ ± 10], 
laminate (Figure 34). The [ ± 30], laminate exhibits more nonlinear response than the 
[ ± 20], laminate because of higher magnitudes of a, (t„) and a greater dependence on the 
transverse and shear response. Due to an even greater dependency upon the transverse and 
shear response, the [ ± 45], laminate shows significant axial softening. 
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4.5.1 .2 Transverse Response 


As can be seen in Figure 35, the agreement between theory and experiment for the trans- 
verse strain response of the four angle-ply laminates is generally good, but not quite as good 
as that for the axial response. In general, the predicted transverse strain response is stiffer 
than the experimental response for all the laminates. The predicted transverse strain re- 
sponse of the [ ± 10], laminate falls between the linear elastic and experiment curves. The 
[ ± 20], laminate's predicted transverse strain response follows the same trend as that of the 
[ ± 10], laminate. The [ ± 30], laminate's experimental transverse strain response shows 
significant softening. Unfortunately, the predicted response does not correlate as well with 
experiment, falling between the linear transverse strain curve and the experimental curve. 
Agreement between theory and experiment is very good for the [ ± 45], laminate in which the 
most softening occurs. Considering the degree of softening, the correlation between the 
predicted transverse strain response and the experimental transverse strain response is ex- 
cellent. 


4.5.2 Response in Material Principal Coordinates 


Examining the predicted stress-strain response in material principal coordinates reveals some 
interesting results. The fiber direction response (Figure 36) is intuitively predictable because 
e\ r is only a function of g 1 (2.11). Therefore, the degree of softening is directly related to the 
magnitude of a v As Figure 36 shows, the magnitude of a, is largest in the [ ± 10], laminate 
and smallest in the [ ± 45], laminate. Therefore, the [ ± 10], laminate exhibits the greatest 
fiber direction softening and the [ ± 45], laminate the least. Unlike the fiber direction re- 
sponse, the stress-strain response in the other principal material directions is not intuitively 
predictable. The non-linear strains in the transverse and shear directions may or may not be 
similar to the softening response in pure transverse loading or shear loading, respectively. 
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The transverse direction response in these laminates can be compared to the transverse di- 
rection response under pure compression loading (see Figure 4 in Chapter 2) and to the linear 
elastic response (Figure 37). As shown in Figure 37, the [ ± 10]„ [ ± 20]„ and [ ± 30], 
laminates show transverse direction stiffening response, and the [ ± 45], laminate exhibits 
softening response. Note that o z is positive in the laminates which show stiffening response 
([ ± 10], , [ ± 20], , and [ ± 30], ), but is negative in the laminate which shows softening 
response ([ ± 45], ). This stiffening response may be due to the endochronic parameters' 
inability to model mixed principal coordinate stress states (a, is negative in all the laminates); 
recall that the parameters were developed from uniaxial tests in which all stresses were of the 
same sign. Another possible explanation of this behavior is that the predicted stiffening re- 
sponse is real and exists within the laminates. However, without experimental results, it is 
difficult to prove or disprove this behavior. As a side note, it is possible that the endochronic 
parameters' inability to model mixed stress states caused the predicted global transverse 
strain response to be stiffer than experiment (Figure 35) for the [ ± 10], , [ ± 20]„ and 
£ ± 30], laminates. Recall that agreement was good for the [ ± 45], laminate in which the 
principal stress states are of the same sign. 

The predicted shear response in these laminates is compared to the pure shear loading re- 
sponse (losipescu test), off-axis loading response, and linear elastic response in Figure 38. 
All predicted curves show a softening response; however, the degree of softening differs 
among the laminates. The [ ± 10], laminate exhibits more shear softening than both the 10° 
off-axis and losipescu tests. At a shear stress of 2.7 ksi, the predicted shear strain deviates 
-7.4% from linear strain, -5.1% from the 10° off-axis test strain, and -8.7% from the losipescu 
test strain. Both the 10° off-axis and losipescu tests show nearly linear response at the low 
shear stress level in the [ ± 10], laminate (less than 3 ksi). The [ ± 20], laminate exhibits 
less shear softening than the 15° off-axis test, but slightly more softening than the losipescu 
test. Actually, within the [ ± 20], laminate, the predicted response is very nearly linear up 
to the maximum level of shear stress (4 ksi). At a shear stress of 3.7 ksi, the predicted shear 
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(a) [ ± 10], Unnotched Laminate 


(b) [ ± 20], Unnotched Laminate 




(c) [ ± 30], Unnotched Laminate (d) [ ± 45], Unnotched Laminate 


Figure 37. or. Versus e, for Unnotched Laminates - AS4/3502 Graphite-Epoxy 
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strain deviates -2.1% from linear strain, 5.5% from the 15° off-axis test strain, and -2.1% from 
the losipescu test strain. The shear response in the [ ± 30], laminate is softer than in the 
30° off-axis test, but slightly stiffer than in the losipescu test. However, the difference between 
the losipescu test shear response and the predicted shear response is very small. At 7.4 ksi, 
the predicted shear strain deviates -11.4% from linear strain, 3.3% from the 30° off-axis test 


strain, and -3.0% from the losipescu test strain. The £ ± 45], laminate shows more shear 


softening response than the losipescu test but not by a large degree. The 45 s off-axis test 


does not load to a high enough stress level to make a useful comparison with predicted re- 
sults. A useful feature of the [ ± 45], laminate is that the shear response can be extracted. 
A plot of the shear response within the tested [ ± 45], laminate is included in Figure 38. 
Comparison between the finite element prediction and the experimental curve reveals excel- 
lent agreement. At 12.5 ksi, the predicted shear strain deviates -36.4% from linear strain, 
-6.5% from the losipescu test strain, and only 4.4% from the Q ± 45], laminate measured 
strain. 


The previous comparison shows how dependent the predicted shear response is on the stress 
component ratios. This phenomenon is also obvious from the off-axis unidirectional tests. 
The ratios of the principal stress components greatly influence the shear stress-strain profile. 
However, as noted previously, unlike the off-axis unidirectional tests, the [ ± 10]„ [ ± 20]„ 
and [ ± 30], laminates exhibit material principal stress components of opposite signs. Be- 
cause the endochronic constants were not derived in this stress environment, the shear re- 
sponse for these laminates is questionable. It is encouraging, however, that within the 
[ ± 45], laminate (which has principal stresses of the same sign), the experimental shear 
response and predicted shear response agree very well. 
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(a) [ ± 10], Unnotched Laminate 


(b) C ± 20], Unnotched Laminate 




(c) [ ± 30], Unnotched Laminate 


(d) [ ± 45], Unnotched Laminate 


Figure 38. a. Versus £, for Unnotched Laminates • AS4/3502 Graphite-Epoxy 
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4.6 Notched Angle-Ply Laminates 


4.6.1 Far-Field Axial Response 


As mentioned previously, the finite element enslysis for Isminstes with holes did not Iced to 


the experimentally observed failure load in all cases. This was true of the [( ± 10) 12 ], and 
[( ± 20) 12 ], laminates with 1/4" diameter holes. However, as Figure 39 shows, far-field axial 
stress-strain correlation between theory and experiment for these laminates and for the 
[( ± 30) 12 ], and [( ± 45) 12 ], laminates is very good. 


There are some interesting differences to note between the far-field axial response of the 
notched and unnotched laminates. Unlike the response of the unnotched laminates 
(Figure 34), the notched [( ± 20) 12 ], laminate (Figure 39) shows more axial softening re- 
sponse than the notched [( ± 10) 12 ], laminate. The cause for this difference is likely due to 
the stress concentration around the hole which causes non-linear strains at a lower load level. 
Unlike the unnotched [ ± 30], laminate (Figure 34), the notched [( ± 30) 12 ], laminate with a 
1/2" (Figure 39) hole shows very little axial softening response. In fact, the predicted results 
show slightly more softening response than the experimental results. The [( ± 30) 12 ], lami- 
nate's small degree of softening response is not expected and the cause is unknown. On the 
other hand, the [( ± 45) 12 ], laminates with 1/2" holes exhibit significant softening response. 
The predicted axial response of the [( ± 45) 12 ], laminate correlates extremely well with ex- 
perimental results. Due to the degree of softening in this laminate, the excellent correlation 
is very encouraging. 
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4.6.2 Axial Strain Response Near The Hole 


This section focuses on the axial strain at the outer surface of the laminates near the hole 
edge. All experiments were conducted with an axial strain gage located as near to the hole 
edge as possible. The location of the strain gage is illustrated in Figure 40. 


In order to compare finite element results with experiment, the far field axial stress is ploiied 


against the axial strain at 0 = 0°, Z = 11H/12 (the closest gauss plane to the free surface), and 


at various distances from the hole edge in Figure 41 through Figure 44. In order to compare 
1/4" hole laminates with 1/2" hole laminates, the distance D from the hole edge is normalized 
with respect to the radius A. 


The finite element stress-strain results agree satisfactorily with experiment between D/A = 
0.17 and D/A = 0.36 for the [( ± 10) 12 ], (Figure 41) and [( ± 20) 12 ] s (Figure 42) laminates. 
The strain calculated at the gauss point nearest the hole (D/A = 0.0048) is much larger than 
the experimentally measured value. At D/A = 0.17440, the predicted response (Figure 41 and 
Figure 42) correlates reasonably well with experiment; however, at this D/A, the stress-strain 
response is softer than the experimental response. The best correlation for both laminates 
appears to occur at 0.17 < D/A < 0.36. 

The predicted results for the [( ± 20) 12 ], laminate (Figure 42) exhibit an interesting 
phenomonen in that the most non-linear stress-strain response occurs not at the gauss point 
closest to the hole edge, but at the next gauss point away (D/A = 0.02). In fact, the response 
at D/A = 0.0352 is also softer than the response at D/A = 0.0048. The [( ± 10) 12 ], laminate 
does not exhibit this behavior. 


The predicted response of the [( ± 30) t2 ], laminate (Figure 43) correlates very well with ex- 
periment at the first three gauss points nearest the hole (D/A = 0.0044, D/A = 0.02, and D/A 


Laminate Response 


127 



Oy FAR FIELD (KSI) 



Figure 41. Far Field Axial Stress Versus Axial Strain on the Hole Edge for [ ± 10 ] Laminate with 
a 1/4" Hole 
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Figure 42. Far Field Axial Stress Versus Axial Strain on the Hole Edge for [ ± 20 1 Laminate with 
a 1/4" Hole 
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Figure 43. Far Field Axial Stress Versus Axial Strain on the Hole Edge for [ ± 30 ] Laminate with 
a 1/2" Hole * 
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Figure 44. Far Field Axial Stress Versus Axial Strain on the Hole Edge for [ ± 45 1 . Laminate with 
a 1/2" Hole 
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= 0.0356). Of these three gauss points, the results at the third gauss point (D/A = 0.03520) 
appear to agree best with experiment. Similar to the predicted response in the [( ± 20) 12 ], 
laminate, the [( ± 30) 12 ], laminate's stress-strain response at D/A = 0.02 is slightly softer 
than the response at D/A = 0.0044, but this behavior is not as pronounced in the [( ± 30) 12 ], 
laminate. In fact, the curves from the nearest gauss point from the hole (D/A = 0.0044) to the 
furthest gauss point from the hole (D/A = 0.1108) all tend to be grouped closer than in the 
[( ± 20) 12 ], and [( ± 10) 12 ], laminates. 

The response of the [( ± 45) 12 ] t laminate (Figure 44) correlates reasonably well with exper- 
iment. Note that the predicted responses at all six gauss points (from D/A = 0.0044 to D/A 
= 0.1108) are nearly identical. The response at D/A = 0.02 and 0.0356 is softer than at D/A 
= 0.0044, but this behavior is nearly indistinguishable. 

It appears from the above results that the laminates with larger layup angles <p show the 
non-linear region adjacent to the hole to extend further towards the free edge. The reason the 
predicted D/A values are larger in the 1/4" hole laminates than in the 1/2" hole laminates may 
be due to their smaller radius value A. At approximately the same distance from the hole 
edge, the D/A value of the 1/4" hole laminates is twice that of the 1/2" hole laminates. 


4.6.3 Stress-Concentration Profile At 9 = 0 ° 


This section examines stress concentration profiles along 0 = 0° from the hole edge to a dis- 
tance D/A = 0.6. Three stress components are examined -- <r e , and a,. At the hole edge, 
due to free surface boundary conditions, these are the only nonzero stresses. Stress profiles 
from linear elastic predictions are compared with those from endochronic theory predictions 
in Figure 45 through Figure 51. The linear results were obtained from the first load step of 
the solution process, while the nonlinear results correspond to the last load step in the loading 
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history. All laminate stresses are normalized with respect to the absolute value of laminate's 
far field stress |cr y | at the corresponding load step. 

4.6.3.1 o 0 Stress Concentration Profiles 

Stress concentration profiles for cr 9 at the three gauss planes closest to the mid-plane 
(Z = H/12, Z = H/4, Z = 5H/12) are shown in Figure 45 through Figure 47 for the four notched 
laminates. At all three gauss planes, the [( ± 10) 12 ], laminate exhibits the highest o e stress 
concentration and the [( ± 45) 12 ], laminate the lowest. Linear elastic analysis predicts a de- 
crease in the peak cr 9 stress concentration as the layup angle <p increases from 10° to 45°. 
The endochronic theory predicts a similar decrease in the peak a 9 stress concentration as the 
layup angle <p increases from 10° to 45° at Z = 5H/12, but at Z = H/4 and Z = H/12, it predicts the 
peak o e stress concentration to be higher in the [( ± 30) 12 ], laminate than in the [( ± 20) 12 ], 
laminate. In the [( ± 10) 12 ] t and [( ± 45) 12 ] s laminates, as intuitively expected, linear elastic 
analysis predicts higher peak a 9 stress concentrations than the endochronic theory. However, 
in the [( ± 20) 12 ], and [( ± 30) 12 ], laminates, linear elastic analysis predicts lower peak a 0 
stress concentrations than the endochronic theory, which is not intuitively expected. At the 
hole edge, the endochronic theory predicts significantly lower o 9 stress concentrations than 
linear elastic analysis at all gauss planes within all of the laminates. Both analyses predict 
the peak o 9 stress concentration to occur in the [( ± 10) 12 ], laminate at D/A = 0 and 
Z = 5H/12. At this point, Figure 47 shows the linear elastic and endochronic theory o 9 stress 
concentrations to be approximately -4.6 and -4.0, respectively. In contrast, the highest peak 
o 9 stress concentration in the [( ± 45) 12 ], laminate occurs at D/A = 0.11120 and Z = 5H/12 and 
has linear elastic and endochronic theory magnitudes of approximately -2.3 and -1.8, respec- 
tively. It is especially interesting to note the shape of the o 9 profiles. With the exception of the 
[( ± 10) 12 ], laminate, at all gauss planes, the maximum stress concentration occurs not at the 
hole edge, but at a small distance inward. This phenomenon is especially pronounced in the 


Laminate Response 


134 









ENDOCHRONIC THEORY LINEAR THEORY 




Figure 47. a e Versus D/A for Notched Laminates with w/d Ratios of 10, Z=5H/12 


Laminate Response 


137 









ENDOCHRONIC THEORY 


LINEAR THEORY 



Laminate Response 


139 







Figure 51. 


O t Versus D/A for Notched Laminates with w/d Ratios of 10, Z = 7H/12 
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endochronic theory results. Finally, note that at D/A = 0.6, all laminates at all three gauss 
planes exhibit approximately the same o 9 stress concentration magnitude of -1.5. At larger 
D/A values, the o, stress concentration approaches the value -1 because a 9 approaches the 
magnitude of o y . 

4. 6. 3. 2 t 0z Stress Concentration Profiles 

The interlaminar t ez shear stress concentration profiles at gauss planes directly below and 
above the ply interface (Z = 5H/12, Z = 7H/12) are shown in Figure 48 and Figure 49. These 
figures show the [( ± 30 ) 12 ], laminate to have the highest t te stress concentration and the 
[( ± 10 ) 12 ], laminate the lowest. In all laminates, linear elastic analysis predicts significantly 
higher peak t 9z stress concentrations than the endochronic theory. This is especially pro- 
nounced in the [( ± 20 ) 12 ], laminate. Other than this difference in peak magnitudes, the linear 
elastic and endochronic theory profiles are nearly identical. Note that the boundary layer in 
which t 9z exists is very small; x 9z is essentially zero for D/A ^ 0.08 for all laminates. 

4.6.3. 3 a Stress Concentration Profiles 

z 

The CT t out-of-plane normal stress concentration profiles, at gauss planes directly below and 
above the ply interface, are shown in Figure 50 and Figure 51. These figures show that linear 
elastic and endochronic theory analysis predict nearly identical profiles at points away from 
the hole edge. However, near the hole edge, the profiles differ significantly, especially in the 
[( ± 20) 12 ], laminate. For the [{ ± 20), 2 ] s laminate, the endochronic theory predicts the a z 
stress concentration to be nearly 0.6 on the hole edge and at Z=7H/12 (Figure 51), compared 
to a linear elastic value of approximately zero. This high magnitude of a t may be the result 
of the zero-slope condition (see "Unloading Considerations" at the beginning of this chapter) 
imposed on another out-of-plane stress component. Note the oscillation about a z = 0 pre- 
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dieted by both analyses. Although the oscillation magnitudes predicted by the endochronic 
theory are larger than those of linear elastic analysis near the hole edge, both analyses show 
a similar pattern. These oscillations may or may not be real. Often times, displacement 
based finite element analysis gives stress oscillations in regions of large strain gradients. 


4.6.4 Material Principal Coordinate Stress Profiles Around The Hole 


This section examines profiles for all six principal stress components, at the hole edge, for 
-90° ^ 0 £ 90°. As in the previous section, the stresses are normalized with the absolute 
value of far field stress |o y |. In order to minimize the number of plots, the in-plane stress 
component profiles are examined at the mid-plane of the first ply (Z = H/4) and the out-of-plane 
stress component profiles are examined at the gauss plane directly below the ply interface 
(Z = 5H/12). These principal stress profiles can be seen in Figure 52 through Figure 57. 

4.6.4.1 Axial Stress oi 

Examination of the a, profiles of the four angle-ply laminates (Figure 52) shows the linear and 
non-linear profiles to be very similar, except in the [( ± 45) 12 ] s laminate. The [( ± 45) 12 ], 
laminate shows a difference between linear elastic and endochronic results, not so much in 
the shape of the curve, but in magnitude. The highest a, stress concentration in the 
[( ± 45) 12 ] t laminate is predicted at the same angle for both analyses, but has a value of -5.0 
for linear results and a value of over -7.0 for non-linear results, a difference of 40%. The other 
laminates show relatively small differences in stress concentration. However, note that in 
both the [( ± 10) 12 ], and [( ± 45) 12 ], laminates, the endochronic theory predicts the largest 
peak stress concentration while in the [( ± 20) 12 ], and [( ± 30) 12 ], laminates, linear elastic 
analysis gives the largest peak stress concentration. As is intuitively expected for both linear 
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(a) [( ± 10) 12 ], Laminate - 1/4" Hole 


(b) [( ± 20) 12 ], Laminate - 1/4" Hole 




-90 -46 0 45 90 


-90 -45 0 46 91 


(c) [( ± 30) 12 ], Laminate - 1/2" Hole 


(d) [( ± 45) 12 ], Laminate - 1/2" Hole 


ENDOCHRONIC THEORY 


LINEAR THEORY 
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Figure 52. a, Versus 6 for Notched Laminates with w/d Ratios of 10, Z=H/4 
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(a) [( ± 10) 12 ], Laminate - 1/4" Hole 


(b) [( ± 20) 12 ], Laminate - 1/4" Hole 
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(c) [( ± 30) 12 ] t Laminate - 1/2" Hole 


(d) [( ± 45) 12 ], Laminate - 1/2" Hole 
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Figure 54. o 3 Versus G for Notched Laminates with w/d Ratios of 10, Z = 5H/12 


Laminate Response 







(a) [( ± 10) 12 ], Laminate - 1/4" Hole 


(b) [( ± 20) 12 ] t Laminate - 1/4" Hole 
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(c) [( ± 30)„], Laminate - 1/2" Hole 


(d) [{ ± 45) 12 ], Laminate - 1/2" Hole 
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Figure 55. o 4 Versus 0 for Notched Laminates with w/d Ratios of 10, Z = 5H/12 
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(a) [( ± 10) 12 ], Laminate - 1/4" Hole 


(b) [( ± 20) 12 ], Laminate - 1/4" Hole 
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(c) [( ± 30) 12 ], Laminate - 1/2" Hole 


(d) [( ± 45) 12 ], Laminate - 1/2" Hole 
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Figure 56. o, Versus 6 for Notched Laminates with w/d Ratios of 10, Z=5H/12 
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(b) [( ± 20) 12 ], Laminate - 1/4" Hole 




(c) [( ± 30) 12 ], Laminate - 1/2" Hole 


(d) [( ± 45) 12 ], Laminate • 1/2" Hole 
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Figure 57. a. Versus 0 for Notched Laminates with w/d Ratios of 10, Z=H/4 
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elastic and endochronic theory results, the peak value of a, occurs at approximately -<p (the 
layup angle of each laminate). 

4.6.4.2 Transverse Stress 02 

Examination of the cr 2 profiles of the four laminates (Figure 53) shows the linear elastic and 
endochronic theory profiles to be very similar in the [( ± 10) 12 ] t and [( ± 20) 12 ] t laminates 
and relatively dissimilar in the [( ± 30) 12 ], and [( ± 45) 12 ], laminates. Note that the dissimi- 
larities tend to occur in a very narrow range of 0 angles. Also, the dissimilarities tend to be 
the most pronounced at the largest values of o 2 . The jaggedness in the linear elastic profiles 
is probably due to local stress oscillations at the gauss points which, as mentioned earlier, 
are common to displacement based finite element analysis in regions of high strain gradients. 
The jaggedness in the endochronic theory profiles is likely due to a combination of the afore- 
mentioned oscillations and to the zero-slope condition forced upon stress components if they 
start to unload in the analysis. Specifying the value of a specific stress component disrupts 
equilibrium and causes the other components to shift to re-establish equilibrium. The 
endochronic theory a 2 stress concentration is generally lower than the linear elastic magni- 
tude in all of the laminates except for the [( ± 30) 12 ], laminate. The [( ± 30) 12 ], laminate's 
profiles are very similar except in a region from 0S —30° to 0 = 20°. In this region, the 
endochronic theory stress concentration increases abruptly, likely caused by the aforemen- 
tioned zero-slope imposition on a certain stress component and the need to re-establish 
equilibrium. 

4.6.4.3 Interlaminar Normal Stress 03 

The linear elastic and endochronic theory stress profiles for the a 3 interlaminar stress 
(Figure 54) are very dissimilar for all the laminates except the [( ± 10) 12 ], laminate. In this 
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laminate, the endochronic theory stress magnitudes are slightly higher than the linear elastic 
magnitudes. Note the [( ± 10) 12 ], laminate's nearly anti-symmetric profile about 0 = 0°. The 
[( ^ 20) 12 ], laminate's linear elastic and endochronic theory profiles agree reasonably well 
from 0 = -90° to 0s -20° and from 0S25 0 to 0 = 90°. However, between 0S -20° and 
0 = 25°, the curves show a significant difference. This difference is mainly caused by a large 
jump in magnitude at 0S -5°. Again, this jump is likely caused by the need to re-establish 
equilibrium due to a zero-slope imposition on another out-of-plane component. The 
[( ^ 30) 12 ], and [( ± 45) 12 ], laminates' endochronic theory predicted a 3 stress profiles exhibit 
radical jumps and show stress magnitudes much higher than the linear elastic profiles. Note 
that in all the laminates, except for the jumps in the endochronic theory profiles, the <y 3 stress 
concentration is relatively low, the absolute value not exceeding approximately 0.1. 

4.6.4.4 Interlaminar Shear Stress 04 

The linear elastic and endochronic theory stress concentration profiles for the o 4 (x 23 ) 
interlaminar shear stress (Figure 55) are similar in all of the laminates except for the 
[( ± 20) 12 ] t laminate where the profiles show a significant difference in magnitude. In this 
laminate, the largest endochronic theory stress concentration is nearly 6 times greater than 
the largest linear stress concentration. The [( ± 10) 12 ], laminate shows the endochronic 
theory positive stress magnitude to be larger and the negative magnitude to be smaller than 
the linear elastic stress magnitude. The [{ ± 30) 12 ], laminate's endochronic theory peak 
stress magnitude is larger than the linear elastic peak magnitude. In contrast, the 
[( ± 45) 12 ], laminate's endochronic theory predicted stress profile exhibits smaller magni- 
tudes than linear elastic analysis. The largest predicted stress concentration occurs in the 
[( ± 45) 12 ], laminate at approximately -20°. Linear elastic analysis predicts a value of nearly 
-0.7 while endochronic theory analysis predicts a value near -0.5. The [{ ± 30) 12 ], laminate 
exhibits magnitudes nearly as high. These large magnitudes of <j 4 stress concentration show 
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the potential for delamination. The magnitudes in the [( ± 20) 12 ], and [( ± 30) 12 ], laminates 
are not nearly as high. 


4.6.4.5 Interlaminar Shear Stress a« 

The linear elastic and endochronic theory a 5 (t 13 ) stress concentration profiles (Figure 56) 
agree reasonably well for all the laminates except the [( ± 20) 12 ], laminate (similar to the a 4 
profiles). This laminate shows odd behavior in that the endochronic theory stress concen- 
tration is opposite in sign to the linear elastic stress concentration, except in the region from 
0= —25° to 0 = 25°. In this region, the profiles are of the same sign, but the endochronic 
theory stress magnitude is much larger than the linear elastic magnitude. In the [( ± 10) 12 ]„ 
[( ± 30) 12 ]„ and [( ± 45) 12 ], laminates, the endochronic theory stress magnitudes are gener- 
ally lower than the linear elastic magnitudes. Again, the largest discrepancies occur near the 
stress concentration peaks. 

4.6.4.6 Inplane Shear Stress oe 

The Oj (t 12 ) linear elastic and endochronic theory stress concentration profiles (Figure 57) are 
similar for all the laminates. In general, they show the endochronic theory stress magnitudes 
to be smaller than the linear elastic magnitudes. Differences are most pronounced at the 
profile peaks. Not only are the linear elastic and endochronic theory stress profiles similar, 
but the profiles for all laminates are very similar in shape, especially in the [( ± 10) 12 ], , 
[( ± 20) 12 ], , and [( ± 30) 12 ], laminates. 
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4.6.5 Cylindrical Coordinate Stress Profiles Around The Hole 


This section examines cylindrical coordinate stress profiles at the hole edge from 
-90° ^ 0 £ 90°. The in-plane normal stress, o e , is plotted at the mid-plane of the first ply 
(Z = H/4) in Figure 58, and the interlaminar shear stress, t*. is plotted at the gauss plane be- 
low the ply interface (Z = 5H/12) in Figure 59. 

4.6.5. 1 Circumferential Stress o 0 


Examination of the o 9 stress profiles (Figure 58) around the hole shows the linear elastic and 
endochronic theory profiles to be very similar except in the [( ± 45) 12 ], laminate. The 
[( ± 45) 12 ], laminate shows the endochronic theory stress magnitude to be significantly higher 
from 0 = 30° to 0S5O°, while from 0= -25° to 0 = 30°, the linear elastic stress magnitude is 
higher. The [( ± 10) 12 ], laminate also shows the endochronic theory peak stress magnitude 
to be higher than the linear peak magnitude by a small amount. However, at all other values 
of 9, the endochronic theory stress magnitude is lower than that predicted by linear theory. 
In both the [( ± 20) 12 ], and [( ± 30) 12 ], laminates, the endochronic theory stress magnitude 
is lower than the linear elastic magnitude for all values of 0. Note how similar the a 9 stress 
profiles are to the o, profiles in the previous section. This most likely occurs because o, is the 
largest stress component magnitude in all of the laminates and is the dominant term in the 
transformation from principal coordinates to cylindrical coordinates. 

4.6.5.2 Interlaminar Shear Stress x Bz 

The t 9i stress profiles (Figure 59) show the maximum value of t 9i to occur at 0sO° for all of 
the laminates. In fact, all of the linear elastic stress profiles have very similar shapes. As <p 
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(a) [( ± 10) 12 ], Laminate ■ 1/4" Hole 


(b) [( ± 20) 12 ], Laminate • 1/4" Hole 




(c) [( ± 30) 12 ], Laminate • 1/2" Hole 


(d) [( ± 45) 12 ], Laminate • 1/2" Hole 
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Figure 58. a Q Versus 0 for Notched Laminates with w/d Ratios of 10, Z=H/4 
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(a) [( ± 10) 12 ], Laminate - 1/4" Hole (b) [( ± 20) 12 ], Laminate - 1/4" Hole 




(c) [( ± 30) 12 ], Laminate - 1/2" Hole 


(d) [( ± 45) 12 ], Laminate - 1/2" Hole 


ENDOCHRONIC THEORY 


--o- LINEAR THEORY 



Figure 59. x 0| Versus 6 for Notched Laminates with w/d Ratios of 10, Z=5H/12 
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(the layup angle of the laminates) increases, the 0 bandwidth of high stresses around the 
hole widens. In fact, the [( ± 45) 12 ], laminate shows a nearly constant peak value of t 0i from 
0s -25°to 0 = 25°. On the other hand, the peak value of t 9i occurs in a very limited region 
around 0 = 0° in the [( ± 10) 12 ], laminate. The endochronic theory stress magnitudes are 
smaller than the linear elastic magnitudes for nearly all values of 0 in all of the laminates. 
Only in the [( ± 10) 12 ], laminate is the endochronic theory stress magnitude larger than the 
linear elastic magnitude from 0 s 15° to 0 s 45°. The reduction in the peak stress concentration 
due to endochronic theory analysis is as great as 30% in the [( ± 45) 12 ], laminate. 


4.7 Summary 


The overall stress-strain correlation between endochronic theory finite element analysis and 
experiment is very good. The endochronic theory was shown to predict axial and transverse 
response more accurately than linear elastic analysis in all of the laminates, especially for the 
[ ± 45], and [( ± 45) 12 ] s laminates. However, the predicted transverse response was not as 
good as the predicted axial response in the [ ± 10]„ [ ± 20], , and [ ± 30], unnotched 
laminates (where <j, and a 2 are of opposite sign). A possible cause for the less accurate 
correlation in these laminates is that the endochronic parameters are unable to accurately 
model mixed stress states. 

In spite of the generally good correlation, the endochronic theory finite element analysis is 
limited by numerical problems. The endochronic constitutive equations predict unloading of 
certain stress components in the angle-ply laminates at some point during their load history. 
This, coupled with the constant stress and constant strain algorithms which will not converge 
near unloading points, limits the utility of the finite element program AMNISAP. The unloading 
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problem is circumvented somewhat by the zero-slope condition. However, this is only a 
"patch-up* solution and results generated after unloading occurs are questionable. 

The comparison between stresses from linear elastic and endochronic theory analyses 
reveals some interesting results. In general, the stress concentrations from endochronic 
theory analysis are smaller than those from linear elastic analysis. However, at various lo- 
cations, the endochronic theory predicts stress magnitudes greater than linear elastic analy- 
sis. The endochronic theory stress profiles are generally more jagged than the linear elastic 
profiles, sometimes jumping to extreme magnitudes. An explanation for this unexpected 
phenomonen is that another stress component at the same location unloaded and the zero- 
slope condition was imposed. To re-establish equilibrium, the observed stress component 
shifted to a higher than expected value. 
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5.0 Failure 


5.1 Introduction 


In a previous analysis of compression loaded angle-ply laminates with central holes, Burns 
et al 1 found ultimate failure load predictions to be very conservative in comparison to exper- 
imental results. However, Burns et al prediction's of the mode and location of failure were 
quite good. In the present study, the main reason for using a nonlinear constitutive relation- 
ship was to asses the effect of nonlinear material response on the prediction of the failure 
load. It was reasoned that linear elastic analysis predicts stresses around the hole to be much 
higher than they are in real laminates. Inclusion of nonlinear material response in the anal- 
ysis should result in more accurate (and lower) stress levels in the region of stress concen- 
tration around the hole. 

Another possible explanation for the conservative predictions of failure using the finite ele- 
ment model is that the idealized model predicts a stress state near the hole which is much 
higher than the true stress state. This is because the finite element model treats the ply 
interface as being discrete ~ the material properties change discontinuously across the 
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interface. Of course, this mathematical discontinuity is not present in actual laminates; 
therefore, the stresses predicted by the model may be significantly higher than those in real 
laminates * 

This chapter presents a failure analysis of both unnotched and notched laminates based upon 
a nonlinear stress analysis and the tensor polynomial failure criterion. The stress results 
generated by AMNISAP were used in a post-processing program which employs the tensor 
polynomial failure criterion to predict the far field load, location (for laminates with holes), and 
mode of failure. The tensor polynomial was chosen as the failure criterion because it ac- 
counts for stress interaction, which may be especially important in the complex, three- 
dimensional stress field around the hole. Results of the failure analysis based upon linear 
and nonlinear material response are compared to experimental results. 

Application of the tensor polynomial (a point-wise failure criterion) at points around the hole 
edge gives a prediction of the initial (first) failure event in the laminate. Often, this does not 
correspond to the ultimate load carrying capacity in real laminates. Damage, in the form of 
micro-cracks, usually forms in the highly stressed region of the hole edge prior to final frac- 
ture. One approach for the application of a point-wise failure criterion to notched laminates 
is to base failure on the stress state at a small distance d 0 from the edge of the hole. The 
distance d 0 from the hole edge has been called a characteristic distance 48 , meaning that for 
a given material it does not change with the layup of the material. This chapter examines the 
tensor polynomial failure predictions at various radii around the hole to assess the correlation 
of theory and experiment. 

The final objective of this chapter is to examine the effect of the w/d (width over hole diameter) 
ratio on failure. Results from linear elastic analysis are compared with those from inelastic 
analysis. 



5.2 The Tensor Polynomial 


The stress tensor polynomial can be written 

Fj<j ( + FijCTjCTj £ 1 (i = 1,2 6) 

where the stress components in material principal coordinates (a,) have been expressed in 
reduced notation. The F, and F,, are second order and fourth order strength tensors, respec- 
tively, determined from tests on unidirectional laminae. Determining the normal stress inter- 
action constants F 12 , F 13 , and F 23 requires difficult to perform biaxial tests. Fortunately, these 
constants have been previously shown to have little influence on the failure predictions. In 
this study, F 12 , F 13 , and F 23 are set to zero. The required components of the strength tensors 
are determined as follows: 


F| " ifc * " * 1 ' 2 ' 3) 

Oi O | 

F|i = ( i = 1 - 2 - 3 - nosum ) 
(CT| ) (<?i ) 

F n = - J |t Y (i = 4,5,6, no sum) 

(<7j ) 

F, = 0 (i = 4,5,6) 

Fjj = 0 (i * j) 
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5.2.1 Strength Parameters for AS4/3502 Graphite-Epoxy 


The AS4/3502 strength data used for failure prediction (Table 9) are the same as those used 
in Burns et ai \ The constants of and op (i = 1,2,3) represent the principal material direction 
normal strengths in compression and tension, respectively. The constants a}"* (i = 4,5,6) rep- 
resent the principal material direction shear strengths. Since the time of Burns' work, im- 
proved test methods have shown the of and of shear strengths to be considerably higher 
than those used by Burns et al. The 1-2 plane [0] losipescu test used to determine 
endochronic constants in Chapter 2 (Figure 7) failed at a shear stress of 13.5 ksi; however, 
failure did not initiate in the shear region and the failure mode was transverse tension, not 
pure shear. Therefore, the true shear strength of AS4/3502 is believed to be higher than 13.5 
ksi. Walrath and Adams 47 observed shear strengths of 15.5 ksi for AS4/3501-6 (a material very 
similar to AS4/3502) using the losipescu shear tests. Considering the above, the 1-2 and 1-3 
plane shear strength (of and of , respectively) used in this study is 14 ksi. Both Burns et al 
value of 9.4 ksi and this new value of 14 ksi are considered in the failure analysis of notched 
and unnotched laminates and comparisons are made. 


5.3 Unnotched Laminates 


Failure analysis of the unnotched laminates presented few difficulties. Because the two- 
dimensional stress state (edge effects were not considered) is constant throughout the in- 
piane and thickness dimensions, ultimate fracture was assumed to occur when failure was 
first predicted by the tensor polynomial. Comparisons of failure based upon linear elastic 
stress analysis and nonlinear stress analysis are shown in Table 10 where the axial failure 
stress and failure strain are given as of and ef, respectively. 
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As can be seen in Table 10, both linear theory and endochronic theory failure predictions for 
the unnotched angle-ply laminates were quite good. Since the tensor polynomial is based on 
stress, the endochronic theory and linear theory predictions of failure stress are nearly iden- 
tical. However, the failure strains predicted by the two theories differ considerably - the 
endochronic theory failure strains correlate with experiment significantly better than the linear 
theory failure strains. 

As shown in Table 10, a shear strength of 9.4 ksi gives conservative failure predictions (both 
of and ey") for all of the laminates, except the [ ± 20], laminate. A shear strength of 14 ksi 
also gives conservative failure predictions except for the [ ± 20], laminate and a;" of the 
[ ± 45], laminate; however, the failure predictions are considerably more accurate (except 
for the [ ± 20], laminate). Using the 14 ksi shear strength, the predicted ay" stresses and 
ey" strains are within 6% of experimental values for the [ ± 10], and [ ± 30], laminates. The 
endochronic theory predicted ay 1 ' and ey H are within 12% and 40% of experiment, respectively, 
for the [ ± 45], laminate. The sy" strain is a significant improvement over that predicted using 
9.4 ksi as the shear strength. The tested [ ± 45], laminate did not fail at the stress and strain 
levels given in Table 10 -- the amount of axial deformation required to cause final fracture of 
this specimen was larger than the travel permitted by the loading frame. The ay" and ey" val- 
ues presented are the maximum values attained during the test. 

The reason for the [ ± 20], laminate's poor correlation between theory and experiment is not 
understood. It is possible that the test specimens contained flaws which caused them to fail 
early. However, this is unlikely because both specimens failed at nearly identical loads. An- 
other possibility is that the [ ± 20], laminate failed due to interlaminar stresses which were 
not modeled in the stress analysis. If this is true, the [ ± 20], laminate's failure mode is ex- 
clusive from the failure modes of the other angle-ply laminates since failure predictions are 
satisfactory for them. 
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Table 10. Comparison of Failure Predictions with Experiment for Unnotched Laminates 


Laminate 

I ±10], 

[±201. 

[±30J. 

[±45], 


G u,t 

Uy 

(ksi) 

n 

■LSI 

■HUB 

m 

E9 

■HI 

mm 

MtiSm 

— 

MSB 

mm 

WEBMl 

of* = of 
= 9.4 ksi 









Linear 

Theory 

-92.658 

-0.5153 

-48.503 

-0.3679 

-27.538 

-0.3725 

-19.058 

-0.6531 

Endochronic 

Theory 

-92.880 

-0.5461 

-48.716 

-0.3784 

-27.363 

-0.3808 

-19.127 

-0.7856 

of = of 

= 14 ksi 









Linear 

Theory 

-95.504 

-0.5310 

-54.183 

-0.4109 

-34.594 

-0.4680 

-30.272 

-1.0374 

Endochronic 

Theory 

-95.031 

-0.5598 

-53.626 

-0.4184 

-33.362 

-0.4759 

-30.205 

-1.9361 

Experiment - 

-96.985 

-0.5683 

-42.490 

-0.3306 

-34.340 

-0.5028 

-26.760 

-1.5910 

Range 

10.630 

(ksi) 

0.0668 

(%) 

m 

0.0232 

(%) 

1 test 

1 test 

1 test 

1 test 


* Note: 

Experimental values are averaged. 

Range indicates the largest difference between experimental values. 
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5.4 Notched Laminates 


Predicting the final fracture load of laminates containing notches is difficult because regions 
near the notch generally fail before the laminate as a whole fails. This local failure may or 
may not be observable. Micro-cracking is usually difficult to observe whereas shear crippling, 
fiber breakage and delamination are much easier to detect. The angle-ply laminates in this 
study displayed no signs of observable (x-ray) damage prior to ultimate failure (with the ex- 
ception of the [( ± 45) 12 ], laminate with a 1/2" hole). Another difficulty in predicting final 
fracture is obtaining an accurate stress state around the hole. As mentioned earlier, the finite 
element mesh models ply interfaces as being discrete. Therefore, the finite element analysis 
may give stresses that are too high. 

Despite these difficulties, predicting the fracture load is only part of understanding the failure 
process of the notched angle-ply laminates. Predicting the mode of failure and the location 
at which failure initiates are important as well. Discussions of fracture, the mode, and the 
location of failure initiation are presented in the following two sections. 


5.4.1 Fracture 

Theory and experiment are compared for laminates with w/d ratios of 10 in Table 11, where 
a“ tt and e“ u represent the far-field axial stress and strain at the fracture load. The fracture load 
is equivalent to the load at which the tensor polynomial predicts failure at the hole edge. As 
mentioned previously, because local failure at the hole edge does not necessarily correspond 
to fracture and because the finite element stresses may be too high, the theoretical o" 1 * and 
e“ tt values in Table 11 are very conservative. A shear strength of 14 ksi gives slightly less 
conservative failure predictions than the 9.4 ksi shear strength; however, the improvement is 


not as pronounced as with the unnotched laminates. The higher value of shear strength im- 
proves the [( ± 45) 12 ], laminates ay" stress most, bringing it to within approximately 50% of 
the experimental ay" stress. 

Similar to the analysis of unnotched laminates, the linear theory and endochronic theory ay" 
stresses in Table 11 are nearly identical for all of the laminates. As mentioned previously, 
because the tensor polynomial is a stress based failure criterion and because the stresses 
are relatively proportional, this is expected. However, unlike the unnotched laminate analysis, 
the linear theory and endochronic theory ey" strains are also nearly identical. The reason the 
ey" strains are so similar is that at the low predicted values of ay", the laminates have not been 
loaded significantly into their softening regions. If the theoretical ay" stresses were close to 
the experimental ay" stresses, the linear theory and endochronic theory ey" strains would differ 
significantly, the endochronic strains being the more accurate. 

From Table 11, it appears that the problems of applying the tensor polynomial at the hole 
edge and finite element stresses that are possibly too high have caused the theoretical ay" 
failure stresses to be very conservative. Because of the low ay" stresses, the endochronic 
theory was not able to show an advantage over linear elastic analysis in predicting failure 
stresses ay" and failure strains ey" . The results from Chapter 4 clearly show that with more 
accurate values of ay" stresses, the endochronic theory would give more accurate predictions 
of ey" strains than linear theory. 


5.4.2 Failure Mode and Location 

The failure mode and location as predicted by linear and endochronic stress analyses agree 
quite well with experiment. Table 12 and Table 13 list the predicted location, fracture stress 
ay" and strain ey", and the tensor polynomial terms at failure for the four angle-ply laminates 
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Table 11. Notched Laminate Comparison of Failure Predictions with Experiment 


Laminate 


mm 

■nsEDm 


—m ii mm 


m 

— 


mm 

■m 

■9 

KB 

■El 

mm 

BK 

KB 

— H 

Wsflm 

a!- 1 ' = CTj" 
= 9.4 ksi 









Linear 

Theory 

-10.696 

-0.0596 

-7.956 

-0.0605 

-7.208 

-0.0973 

-7.952 

-0.2722 

Endochronic 

Theory 

-10.640 

-0.0592 

-7.969 

-0.0606 

-7.222 

-0.0975 

-7.871 

-0.2688 

CT U« - a ul« 

= 14 ksi 









Linear 

Theory 

-14.603 

-0.0812 

-11.072 

-0.0841 

-9.570 

-0.1292 

-11.437 

-0.3906 

! 

Endochronic 

Theory 

-14.929 

-0.0833 

-11.223 

-0.0855 

-9.732 

-0.1315 

-12.195 

-0.4239 

Experiment* 

-76.255 

-0.4054 

-45.010 

-0.3550 

-26.323 

-0.3521 

-25.227 

-1.4387 

Range 

1.350 

(ksi) 

0.0161 

(%) 

0.420 

(ksi) 

0.0057 

(%) 

1.510 

(Ksi) 

0.0289 

(%) 

4.28 

(ksi) 

0.3774 

(%) 


* Note: 

Experimental values are averaged. 

Range indicates the largest difference between experimental values. 
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using the two material models. Figure 60 through Figure 63 show the tensor polynomial 
profile around the hole edge for both linear and endochronic stress analyses. Note that the 
endochronic theory gives a profile nearly identical to the linear elastic profile. In the sections 
that follow, theory and experiment are compared. 

5.4.2.1 10 ° Angle-Ply Laminate With A 114 " Hole 

The [{ ± 10 ) 12 ] b laminate failed in two modes, one a compressive failure across the horizontal 
axis of symmetry, and the other a fracture at 10 0 from the loading axis \ Both of these failures 
initiated at 0 = 0°. As shown in Table 12, both linear and endochronic theory analysis predict 
failure to occur at the first gauss plane above the ply interface (Z = 7H/12) and at an angle of 
0.0 ° from the x-axis. The largest term in the tensor polynomial is the out-of-plane shear term 
F55 T i3» which has a value of 0.51 for both types of analysis. This term, along with the relatively 
high in-plane shear term of F 68 t^ 2 = 0.21 likely represents the shear crippling type failure ob- 
served in the test specimen. The compressive failure mode is harder to explain. The F 1 a 1 and 
terms are relatively small; however, the value of their sum is significant, being approxi- 
mately 0.26. Both terms represent compressive failure, so their sum is the important value. 
Although this value is not as large as the F 55 t^ 3 term, it is the next largest value and thus re- 
presents compressive failure as a second possible failure mode. Note that the linear elastic 
and endochronic theory tensor polynomial terms are nearly identical. Apparently, the stress 
states predicted at failure are very similar for both analyses. Using a shear strength of 14 ksi, 
Table 12 shows the endochronic prediction for failure location to be the same as the elastic 
prediction; however, the tensor polynomial terms are different from those using using a shear 
strength of 9.4 ksi. The compressive failure mode terms (F^ and are higher, while the 
shear mode terms (F 55 t^ and F 66 xy are lower. In fact, the sum of F^ and F^af is approxi- 
mately equal to F 55 t? 3 , both values being roughly 0.40. This is consistent with both failure 
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Table 12. Tensor Polynomial Terms for Laminates with 1/4" Holes 



l(±10)„]. 

U±20)J. 

of = erf 

Linear 

Endochronic 

Linear 

Endochronic 

= 9.4 ksi 

Theory 

Theory 

Theory 

Theory 

Z 

7H/12 

7H/12 

7H/12 

7H/12 

e 

0.0° 

0.0 0 

-2.5 0 

-2.5 0 

WmStmSsMm 

-10,695 

-10,640 

-7,956 

-7,969 

M1IW 

-.0596 

-.0592 

-.0605 

-.0606 

FiCTi 

0.1583 

0.1579 

0.0801 

0.0799 


0.1005 

0.1000 

0.0257 

0.0256 

F 2 a 2 

0.0008 

0.0006 

0.0362 

0.0358 

F 22 cy 2 

0.0000 

0.0000 

0.0004 

0.0004 

^3^3 

-0.0096 

-0.0098 

0.0178 

0.0173 

F 33°3 

0.0000 

0.0000 

0.0001 

0.0001 

wmm 

0.0195 

0.0197 

0.0822 

0.0819 

mmm 

0.5174 

0.5187 

0.5582 

0.5525 

B9I 

0.2130 

0.2111 

0.1993 

0.1965 

a ult = a uH 

Linear 

Endochronic 

Linear 

Endochronic 

= 14 ksi 

Theory 

Theory 

Theory 

Theory 

Z 

7H/12 

7H/12 

7H/12 

7H/12 

e 

0.0 0 

0.0° 

-2.5 0 

-2.5 0 | 


-14,603 

-14,929 

-11,072 

-11,223 

103^91 

-.0812 

-.0833 

-.0841 

-.0855 

F,a, 

0.2173 

0.2186 

0.1113 

0.1119 


0.1893 

0.1917 

0.0497 

0.0502 

F 2 a 2 

0.0012 

0.0015 

0.0503 

0.0532 

F 22 a| 

0.0000 

0.0000 

0.0008 

0.0009 

F 3 cr 3 

-0.0131 

-0.0129 

0.0247 

0.0279 

F 33 a^ 

0.0001 

0.0001 

0.0002 

0.0002 

EH 

0.0368 

0.0380 

0.1587 

0.1614 


0.4027 

0.3993 

0.4453 

0.4384 

EH 

0.1658 

0.1604 

0.1590 

0.1529 
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Table 13. Tensor Polynomial Terms for Laminates with 1/2” Holes 


a u!t = c un Linear 

= 9.4 ksi Theory 


7H/12 


-5.0° 


-7,208 


-.0973 


0.0634 


0.0161 


0.0552 


0.0010 


0.0515 


0.0008 


0.1494 


0.4154 


0.2472 


a ui« = a ui« Linear 

= 14 ksi Theory 


7H/12 


-5.0 0 


-9,570 


-.1292 


0.0840 


0.0283 


0.0731 


0.0017 


0.0683 


0.0015 


0.2623 


0.3015 


0.1794 


Endochronic 

Theory 


7H/12 


-5.0 0 


-7,222 


-.0975 


0.0631 


0.0160 


0.0557 


0.0010 


0.0516 


0.0008 


0.1487 


0.4114 


0.2429 


Endochronic 

Theory 


5H/12 


-2.5 


-9,732 


-.1315 


0.0663 


0.0176 


0.0516 


0.0008 


0.0608 


0.0012 


0.3954 


0.2470 


0.1637 


K±45) 12 ]. 

Linear 

Endochronic 

Theory 

Theory 

5H/12 

5H/12 

-17.5 0 

-17.5 0 

-7,952 

-7,871 

-.2722 

-.2688 

0.0743 

0.0741 

0.0221 

0.0220 

-0.2451 

-0.2414 

0.0188 

0.0182 

0.0534 

0.0570 

0.0009 

0.0010 

0.0906 

0.0898 

0.2648 

0.2670 

0.7202 

0.7135 

Linear 

Endochronic 

Theory 

Theory 

7H/12 

7H/12 

2.5 0 

15.0 0 

-11,437 

-12,195 

-.3906 

-.4239 

0.0505 

0.0265 

0.0102 

0.0028 

-0.5591 

-0.6354 


0.0977 


0.0779 


0.0019 


0.6390 


0.1295 


0.5524 


0.0785 


0.0019 


0.9066 


0.0551 


0.4297 
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TENSOR POLYNOMIAL 
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Tensor Polynomial Profile Around The Hole Edge - [(±10),,]. Laminate With 1/4" Hole, 
Z = 7H/12 
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Figure 61. Tensor Polynomial Profile Around The Hole Edge - [(±20) 12 ] Laminate With 1/4" Hole, 
Z = 7H/12 
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Figure 82. Tensor Polynomial Profile Around The Hole Edge - [(±30),,). Laminate With 1/2" Hole, 
Z = 7H/12 
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TENSOR POLYNOMIAL 
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Figure 63. Tensor Polynomial Profile Around The Hole Edge - [(±45)„). Laminate With 1/2* Hole, 
Z = 5H/12 
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modes being present in the test specimen. This result also lends support for using the higher 
value of shear strength. 


Figure 60 shows that not only are the linear elastic and endochronic theory stress states 
nearly identical at the point of failure, but around the profile of the hole as well. The 
endochronic theory results, represented by a dashed line, are nearly indistinguishable from 
the linear elastic results. The shape of the tensor polynomial profile is very interesting. The 
majority of the profile lies under the value of 0.5. At 9 — 0°, there is a very sharp peak where 
the tensor polynomial reaches 1.0. Note also that the profile is nearly symmetric about 
0 = 0°. This is because of the relatively small fiber angle of 10 °. As the other figures show, 
the near symmetry about 0 = 0° is less pronounced for larger fiber angles. 

5 .4.2.2 20 0 Angle-Ply Laminate With A 114" Hole 

The [( ± 20) 12 ], laminate failed due to a fracture at 20 ° from the loading axis \ Like the 
[( ± 10 ) 12 ]. laminate, this fracture was parallel to the fiber direction in half of the plies and 
fractured fibers in the alternate plies. Failure initiated at 0 = 0°. Table 12 shows that both 
linear elastic and endochronic theory analysis predict failure to initiate in the first gauss plane 
above the ply interface (Z = 7H/12) and at 0 = —2.5° . The dominant tensor polynomial term 
is again the out-of-plane shear term F 55 Tf 3 with a value of approximately 0.55. The next largest 
term is the in-plane shear component F 66 t ? 2 with a value of approximately 0.20. Together, 
these two terms have a sum of 0.75, which represents a shear failure as being the dominant 
failure mode. Using a shear strength of 14 ksi, Table 12 shows a slight redistribution in the 
values of the tensor polynomial terms that lessens the dominant effect of the F„t ? 3 and F 66 t ? 2 
terms and increases the values of the F 1 ct 1 , F 13 o 3 , and F^ifj terms. However, the F 55 Tf 3 term is 
still the largest term at a value of approximately 0.44. As with the [( ± 10) 12 ], laminate anal- 
ysis, the tensor polynomial terms are nearly identical for both linear and endochronic theory 
analysis. 
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Figure 61 shows the tensor polynomial profile around the hole of the [( ± 20) 12 ], laminate to 
be nearly identical to the profile of the [( ± 10) 12 ], laminate. The peak value occurs at 
6 = -2.5° and the peak is not quite as sharp as in the [( ± 10) i2 ], laminate. Also, the near 
symmetry noted in Figure 60 is somewhat more skewed in Figure 61. 

5 . 4 . 2.3 30 ° Angle-Ply Laminate With A 112" Hole 

The [( ± 30) 12 ], laminate failed due to a fracture at 30 0 from the loading axis \ As in the 
[( ± 10) 12 ], and [( ± 20 ) 12 ] s laminates, this fracture was parallel to the fiber direction in half 
of the plies and fractured fibers in the alternate plies. Failure appears to have initiated at 
0 = 0°, although the failure band was fairly wide, thus it is difficult to pinpoint the exact location 
at which failure begins (see Figure 2.25 in Burns et al ’J. Using a shear strength of 9.4 ksi, 
Table 12 shows that both linear theory and endochronic theory analyses predict failure to in- 
itiate in the gauss plane above the ply interface and at 0 = —5.0°. This location is reasonably 
close to that observed in the test specimen. Like the [( ± 20) 12 ], laminate, the dominant 
tensor polynomial term is F 55 Tf 3 with a value of 0.41. The second and third largest terms are 
FhjT^ and F^ifj, respectively. The sum of these shear terms is greater than 0.8, which shows 
that shear failure should be dominant. Note that all of the tensor polynomial terms agree quite 
well between linear elastic and endochronic theory analysis. Using a shear strength of 14 
ksi, Table 13 shows a difference in the location of failure for the two types of analysis. Linear 
analysis predicts failure to occur at the same point <Z = 7H/12, 0=-5.O °), whereas the 
endochronic theory predicts failure to occur at the gauss plane directly below the ply interface 
(Z = 51-1/12) at an angle of 0 = -2.5°. The endochronic theory predicts the failure location 
more accurately than linear elastic analysis when the shear strength is set to 14 ksi. Note that 
the tensor polynomial terms differ fairly significantly between the two analyses and have 
shifted from the values using a shear strength of 9.4 ksi. Linear analysis still shows the 
F 55 t? 3 term to be dominant, but to a lessor degree. However, the endochronic theory shows 
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the F 44 T 23 term to be dominant. In both analyses, the sum of the tensor polynomial shear terms 
using shear strengths of 14 ksi and 9.4 ksi are nearly the same; however, the individual terms 
are redistributed. 

Figure 62 shows the tensor polynomial profile around the hole to be very similar to the pro- 
files of the [( ± 10) 12 ], and [( ± 20) 12 ], laminates shown in Figure 60 and Figure 61, respec- 
tively. The peak has a much smoother profile and the bandwidth of points above 0.5 has 
increased. The near symmetry shape is still present, but more skewed. 

5. 4. 2. 4 45 0 Angle-Ply Laminate With A 1/2 ” Hole 

Like the [ ± 45], unnotched laminate, the [( ± 45) 12 ], laminates with 1/2" holes were not 
tested to failure due to the limited travel of the testing machine. However, these laminates 
did show damage in the form of bands of matrix damage, visible under ultrasonic C-scan, 
which formed at ± 45° to the loading axis \ It is nearly impossible to observe the precise 
location of failure initiation from the C-scan, although the damage bands seem to initiate at 
0= ± 15°. Using a shear strength of 9.4 ksi, Table 13 shows both linear elastic and 
endochronic theory analysis to predict the failure location at the gauss plane directly below 
the ply interface (Z = 5H/1 2) and at 0 = 17.5°, which agrees reasonably with experiment. The 
out-of-plane shear terms (F^, and F 55 t 5 3 ) are much less dominant in this laminate. Instead, 
the in-plane shear term F 66 t5 2 dominants with a value of approximately 0.72. Part of the reason 
this value is so high is that the F 2 a 2 term is negative, which allows the other terms to be larger. 
Shear failure is the obvious failure mechanism predicted by both analyses. Again, the tensor 
polynomial terms are nearly identical for both analyses. Using a shear strength of 14 ksi, 
Table 13 shows a large discrepancy in the predicted location of failure between linear elastic 
and endochronic theory analysis. Both analyses predict failure to occur in the gauss plane 
directly above the ply interface (Z=7H/12). However, linear analysis predicts the angle of 
failure initiation to be 2.5° whereas endochronic theory analysis predicts 0 = 15.0°. The 
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endochronic analysis predicts the correct location more accurately than linear elastic analy- 
sis. Also, the tensor polynomial terms differ significantly between the two analyses. In both 
analyses, the dominant term has shifted from F 66 Ti 2 to However, this shift is most pro- 

nounced for the endochronic theory. Note also that the F 2 o 2 term using a shear strength of 
14 ksi is almost three times larger than the F term using a shear strength of 9.4 ksi. 

The shape of the [( ± 45) 12 ] s laminate's tensor polynomial profile is, again, very similar to the 
previously discussed laminates, as can be seen in Figure 63. The profile is anti-symmetric 
to the previous laminate's profiles because failure is predicted at Z = 5H/12 instead of 
Z = 7H/12. The bandwidth of points with magnitudes above 0.5 has increased, and the near 
symmetrical shape is much more skewed. As with the previously discussed laminates, the 
endochronic theory profile is nearly indistinguishable from the linear elastic profile. 


5.4.3 Failure Away From The Hole Edge 


Whitney and Nuismer 43 applied the characteristic distance theory by evaluating a given failure 
criterion at a set distance d„ from the hole edge along the horizontal centerline of the speci- 
men (0 = 0° here). Burns 1 also used this approach in his work and found the distance d„ to 
be dependent on both hole diameter and laminate ply orientation. In other words, the distance 
d 0 was not found to be "characteristic". 

Instead of examining the state of stress at various distances from the hole along 0 = 0°, this 
study applies the tensor polynomial at various radii (or distances from the hole edge) for 
— 90° ^ 0 ^ 90° . Figure 64 and Figure 65 show the fracture stress ay 1 * as a function of 
distance from the hole edge for the four notched laminates with a w/d ratio of 10. Figure 64 
shows the profiles from linear elastic analysis and Figure 65 shows the endochronic theory 
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profiles. Table 14 lists the failure locations for the same four laminates as a function of dis- 
tance from the hole edge. 

Figure 64 and Figure 65 show the failure profiles from linear elastic analysis and endochronic 
theory analysis to be nearly the same. Note that the [( ± 30) 12 ] s and [( ± 45) 12 ], laminates 
both have a local maximum at D/A = 0.0045. In other words, the predicted failure stress is 
lower at D/A = 0.02 than at 0.0045. In fact, for the [( ± 45) 12 ] s laminate, this local maximum 
is the absolute maximum for the data shown. The [( ± 45) 12 ] s laminate's profile appears to 
have nearly a 0 0 slope after the point D/A = 0.02. These figures show that very large D/A 
values are required to predict the experimentally measured strengths. 

As D/A increases, the location of failure also changes as shown in Table 14. Both linear 
elastic and endochronic theory analysis predict the same failure location in all of the lami- 
nates except for the [( ± 45) 12 ], laminate. However, the discrepancy in this case is small and 
only occurs at D/A = 0.02 and D/A = 0.0355. Note that the through-the-thickness Z location 
varies widely, especially for the [( ± 10) 12 ] s and [( ± 20) 12 ], laminates. The 0 locations tend 
to vary less than the Z locations with the exception of the [( ± 10) 12 ] s laminate, where it jumps 
to nearly 90° at D/A = 0.08. In examining this table, keep in mind that the through-the- 
thickness location Z and the angle 0 are not entirely separate. A point located at a given 0 
and Z = 7H/12 has nearly the same stress state as a point located at -0 and Z = 5H/12. 
For example, in the [( ± 30) 12 ], laminate, failure is predicted to occur at D/A = 0.0045, Z = 
7H/12 and 0 = -3.75°. At D/A = 0.02, failure is predicted at Z = 5H/12 and 0 = 3.75° . Al- 
though these points are different, the failure mechanism and mode are likely the same. 
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Figure 64. Predicted Failure Stress Versus Distance Away From Hole Edge Using Linear Analysis 
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Figure 65. Predicted Failure Streas Veraua Distance Away From Hole Edge Ualng Endochronic 
Theory Analysis 
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Table 14. Location of Failure at Various Distances from the Hole Edge 


l(±10) 12 ]. 


t(±20) 12 ]. 


\ 



Linear 

Theory 

Endochronic 

Theory 

Linear 

Theory 

Endochronic 

Theory 


Z 

0 

Z 

0 

Z 

0 

Z 

0 

0.00000 

7H/12 

0.0 

7H/12 

0.0 

7H/12 

-2.5 

7H/12 

-2.5 

0.00450 

5H/12 

-0.85 

5H/12 

-0.85 

7H/12 

-3.75 

7H/12 

-3.75 

0.020000 

H/12 

-3.75 

H/12 

-3.75 

5H/12 

0.85 

5H/12 

0.85 

0.03550 

H/12 1 

-8.35 

H/12 

-8.35 

5H/12 

0.85 

5H/12 

0.85 

0.04900 

7H/12 

-18.75 

7H/12 

-18.75 

5H/12 

0.85 

5H/12 

0.85 

0.08000 

H/12 

89.15 

H/12 

89.15 

5H/12 

0.85 

5H/12 

0.85 

0.11100 

H/12 

89.15 

H/12 

89.15 

11H/12 

3.75 

1 1 H/12 

3.75 


U±30)J. 

[<±45) 12 ], 


Linear 

Endochronic 

Linear 

Endochronic 


Theory 

Theory 

Theory 

Theory 

D/A 

z 

0 

Z 

0 

Z 

0 

Z 

0 

0.00000 

7H/12 

-5.0 

7H/12 

-5.0 

5H/12 

17.5 

5H/12 

17.5 

0.00450 

7H/12 

-3.75 

7H/12 

-3.75 

5H/12 

15.85 

5H/12 

15.85 

0.020000 

5H/12 

3.75 

5H/12 

3.75 

5H/12 

14.15 

5H/12 

8.35 

0.03550 

5H/12 

0.85 

5H/12 

0.85 

5H/12 

8.35 

5H/12 

3.75 

0.04900 

5H/12 

0.85 

5H/12 

0.85 

5H/12 

3.75 

5H/12 

3.75 

0.08000 

5H/12 

3.75 

5H/12 

3.75 

5H/12 

3.75 

5H/12 

3.75 

0.11100 

11H/12 

-0.85 

11 H/12 

-0.85 

1 1 H/12 

-3.75 

11 H/12 

-3.75 
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5.4.4 The Effect Of W/D On Failure 


Thus far, this chapter has only examined failure in those laminates with w/d ratios of 10. 
However, [( ± 10) 12 ] t and [( ± 20) 12 ], laminates with 1/2' holes and w/d ratios of 5 were also 
tested \ This section examines the effect this reduction in w/d has on both the experimental 
and predicted fracture stress, o“", and strain, eJ“. 

If a hole in a plate caused no stress concentration, the reduction in strength due to a decrease 
in the ratio w/d would simply be equal to the reduction in cross-sectional area. However, in 
this study, the angle-ply laminates have in-plane and out-of-plane stress concentrations at 
both the hole edge and at the edge of the laminate. For the anisotropic plates of this study, 
the stress concentrations at the hole edge are dependent on the hole size as well. With this 
in mind, it is not surprising that the reduction in strength does not follow the reduction in 
cross-sectional area. Comparing the laminates with w/d ratios of 10 to those with w/d ratios 
of 5, there is only a 11.1% reduction in cross-sectional area. However, as Table 15 shows, 
experiments reveal a reduction in strength of approximately 18.5% and a reduction in ultimate 
failure strain of approximately 20% for both laminates. The finite element analyses do not 
predict the same strength reduction as experiment. Both linear theory and endochronic the- 
ory analysis predict approximately an 11% reduction in strength and ultimate failure strain for 
the [( ± 10) 12 ], laminates. For the [( ± 20) 12 ], laminates, linear elastic analysis predicts ap- 
proximately a 9% reduction in strength and ultimate failure strain while the endochronic the- 
ory predicts approximately a 10% reduction. Therefore, the finite element results tend to 
follow the reduction in cross-sectional area of approximately 11%. Because the finite element 
results are extremely conservative in predicting the fracture load, it may not be realistic to 
expect them to agree with the experimental strength reduction. 
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Table 15. The Effect of w/d on Notched Laminate Failure 


[( ± 10), 2 ], 

ay'* (ksi) 

£“"(%) 

Z 

6 

Experiment 

w/d = 10 

-76.255 

-0.44485 


0.0 

(Average) 







w/d = 5 

-62.020 

-0.35170 


0.0 

Linear 

w/d = 10 

-10.696 

-0.05960 

7H/12 

0.0 

Elastic 







w/d = 5 

-9.568 

-0.05320 

7H/12 

0.0 

Endochronic 

w/d = 10 

-10.640 

-0.05920 

7H/12 

0.0 

Theory 







w/d = 5 

-9.441 

-0.05240 

7H/12 

0.0 

[( ± 20) 12 ], 

a" 1 ’ (ksi) 

£“"(%) 

Z 

0 

Experiment 

w/d = 10 

-45.010 

-0.35495 


0.0 

(Average) 







w/d = 5 

-36.675 

-0.28365 

— 

0.0 

Linear 

w/d = 10 

-7.956 

-0.06050 

7H/12 

-2.5 

Elastic 







w/d = 5 

-7.226 

-0.05490 

7H/12 

-2.5 

Endochronic 

w/d = 10 

-7.969 

-0.06060 

7H/12 

-2.5 

Theory 







w/d = 5 

-7.190 

-0.05460 

7H/12 

-2.5 
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5.5 Summary 


Failure predictions for the unnotched angle-ply laminates were very good whereas failure 
predictions for the notched laminates were very conservative. The predicted failure response 
of unnotched laminates correlated quite well with experiment using a shear strength of 9.4 ksi. 
However, a shear strength of 14 ksi yielded even better correlation for all laminates except the 
[( ± 20) 12 ], laminate whose failure prediction was non-conservative with the shear strength 
of 9.4 ksi. The endochronic theory predicted ct“ k stresses almost identical to those of linear 
analysis. However, the endochronic theory predicted strains more accurately than linear 
analysis, especially for the [( ± 45) 12 ], laminate which exhibits the most ductility. 

Using the tensor polynomial at the hole edge to predict the fracture load of notched laminates 
gave extremely conservative results. One reason for this is likely due to local failure in the 
region around the hole before ultimate failure occurs. Another is that the predicted stresses 
around the hole are too high due to the finite element model idealization of discrete ply 
interfaces. 

The tensor polynomial predicted the correct location and mode of failure initiation for the 
notched laminates. However, changing the shear strength from 9.4 ksi to 14 ksi showed a 
difference in the predicted location of failure for some of the laminates. Almost all failures 
were predicted to be caused by t 13 or t 12 shear stresses or both of them combined. A shear 
failure mode agrees well with the observed failure mode of the test specimens. 

Endochronic theory analysis gives a stress state nearly identical to that of linear analysis, for 
a given loading, in both the unnotched and notched laminates. Therefore, both analyses 
predicted of* stresses to be nearly the same for all notched laminates. Unlike the unnotched 
laminates, the linear theory and endochronic theory analyses predicted nearly identical 


values as well, the reason being that none of the laminates were deformed significantly into 
their softening regions before failure was predicted. 

In conclusion, the endochronic theory failure analysis of notched laminates correlated with 
experiment only slightly better than the linear analysis due to the limits of the point-wise fail- 
ure criterion and the finite element model that were used. Chapter 4 showed nonlinear ma- 
terial behavior to be an important factor in analyzing notched angle-ply laminates. Therefore, 
without the limitations imposed by the point-wise failure criterion and the discrete ply interface 
finite element model, the endochronic theory would likely predict significantly more accurate 
failure stresses, ct“" , and failure strains, e"'*. 
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6.0 Conclusions 


This study has analyzed the loading response and failure of unnotched and notched angle-ply 
laminates subjected to uniaxial compression. The analysis was performed using the three- 
dimensional material nonlinear finite element program AMNISAP, which utilizes the 
endochronic theory as a constitutive relationship. Analytical results were compared to ex- 
periment and, where appropriate, results from the endochronic theory were compared to re- 
sults from linear theory. The conclusions of this study from the previous chapters are 
presented below. 


6.1 The Endochronic Theory 


The endochronic theory is a macroscopic constitutive theory that was used in this study to 
model AS4/3502 graphite-epoxy. The theory, originally developed by Valanis 11 12 and modified 
to model transversly isotropic materials by Pindera and Herakovich 19 20 , is based on irre- 
versible thermodynamics and uses the concept of internal variables. The transversely 
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isotropic theory was extended to three-dimensions in this study to facilitate the analysis of 
laminates with cutouts. 

Like the 2-D equations, the 3-D endochronic equations are based on two fundamental laws of 
thermodynamics: 1) Conservation of Energy, and 2) The Dissipation (Clausius-Duhem) Ine- 
quality. The 3-D equations contain 18 constants (15 independent) which are unique to 
AS4/3502 graphite-epoxy. These constants can be separated into reversible constants and 
irreversible constants. The determination of the constants is described as follows: 

1. The six reversible constants (five independent) were used in the expression for reversible 
strain and consist of the compliance matrix Ajj. Because the compliance matrix models 
the linear response of the material, it was determined from the initial response of various 
unidirectional lamina tests. Straight lines were least squares fit through the initial data 
of the lamina tests and the slopes of these lines were used to determine the Ay 

2. Twelve irreversible constants (10 independent) were used in the expression for irrevers- 
ible strain; they are the Sy matrix and the n„ , n 2 , B^, n 6 , and Bjj, parameters. These 
constants, which model the nonlinear response of the material, were also determined 
from various unidirectional lamina tests. The logarithm of a strain component was plotted 
against the logarithm of its corresponding stress component and a straight line was least 
squares fit through this data. The slope and strain axis intercept of these lines were used 
to determine the constants. Sixteen sets of the irreversible constants were determined; 
however, only five sets were considered in the analysis - those whose Sy matrix was 
positive-definite. These five sets were used to model off-axis lamina response and 
unnotched angle-ply laminate response. The parameter set which best represented the 
experimental lamina and laminate tests was chosen for use in this study. 
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6.2 The Finite Eiement Program AMNISAP 


The three-dimensional material nonlinear finite element program AMNISAP utilizes the 
endochronic constitutive theory and constant stress and constant strain solution algorithms 
to solve the nonlinear system of equations. It also utilizes a symmetry relationship called in- 
version symmetry. In this study, the program was vectorized to run on NASA's Cyber 205 
super computer. The important features of AMNISAP are outlined below. 


1. The inversion symmetry utilized in AMNISAP consists of rotating a plate through 180° 
about the axis of symmetry. Defining an x, coordinate system (i = 1,2,3), with the x 3 coor- 
dinate the axis of symmetry, a symmetry transformation changes the coordinates 
x a (a = 1,2) of each material point of the plate into — x a . The inversion symmetry analysis 
in this study resulted in considerable storage and run time savings over conventional half 
plate analysis. 


2. The vectorizing utility VAST was used to vectorize the three subroutines COLSOL, ELKAY, 
and STRESS in AMNISAP because they were found to use the majority of CPU time. 
Based on a 20 element mesh, a total speed factor of 12 was realized by using vectorized 
inversion symmetry analysis over non-vectorized half plate analysis. Considerably higher 
speed factors were likely realized for the analysis of notched angle-ply laminates in this 
study. 

3. The constant stress and constant strain solution algorithms used in AMNISAP are based 
on the concept of initial strain and are classified as step-iterative methods. They each 
utilize a psuedo load vector which contains all the nonlinear terms of the equations. 
Therefore, the global stiffness matrix is formed only once and remains constant through- 
out the analysis. Both algorithms were included in AMNISAP because one of them might 
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be better suited for a given problem. The user can choose the algorithm which best fits 
the problem at hand. 


6.3 Laminate Response 


The response of unnotched and notched angle-ply laminates subjected to uniaxial 
compressive loading was analyzed and compared to experiment. All theoretical results were 
obtained using AMNISAP. The highlights of this analysis are outlined below. 


6.3.1 Unnotched Laminates 

1. The unnotched angle-ply laminates' axial and transverse response agreed very well with 
experiment. The agreement of theory and experiment for the [ ± 45], laminate was es- 
pecially good, considering this laminate's substantial softening response. Although the 
theoretical/experimental correlation for transverse response was good, it was not quite 
as good as that for axial response in the [ ± 10], , [ ± 20]„ and [ ± 30], laminates. A 
possible cause for the less accurate transverse response is that the endochronic 
equations used cannot accurately model mixed stress states (a, and o 2 are of opposite 
sign in the [ ± 10]„ [ ± 20]„ and [ ± 30], laminates). 
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6.3.2 Notched Laminates 


1. The endochronic equations were found to predict unloading of various stress components 
in the analysis of the notched laminates. At a critical load level, this unloading phenom- 
enon prevented higher loading by causing the solution to diverge. A condition referred 
to in this study as the zero-slope condition was implemented in AMNISAP to allow loading 
past the critical load level. The zero-slope condition forces an unloading stress compo- 
nent's stress-strain modulus to zero for the remainder of the loading history. 

2. Theoretical/experimental correlation of the notched angle-ply laminates' far field axial 

response was excellent. Due to the unloading phenomenon mentioned above, the 
[( 10), 2 ], and [( ± 20), 2 ], laminates did not reach the experimentally observed failure 

load. 

3. The theoretically predicted axial strain response near the hole agreed well with exper- 
iment. However, the predicted response was dependent on the distance from the hole, 
especially in the [( ± 10), 2 ], and [( ± 20), 2 ], laminates. The D/A value (nondimensional 
distance from the hole) at which the [( ± 10) 12 ], and [( ± 20) 12 ], laminates best agreed 
with experiment was different from the D/A value for the [( ± 30), 2 ], and [( ± 45), 2 ], 
laminates. However, the [{ ± 10), 2 ], and [( ± 20) 12 ], laminates had 1/4" holes whereas 
the [( ± 30), 2 ], and [( ± 45) u ], laminates had 1/2" holes, which may have influenced the 
placement of the strain gage on the test specimens. 

4. The predicted a 9 , t 9l , and o z stress concentration profiles were examined along 0 = 0° 
from the hole edge to a distance D/A =0.6. The endochronic theory profiles differed sig- 
nificantly from the linear theory profiles, especially for the out-of-plane shear stress t 9i . 
The important features of the stress concentration profiles are given below: 
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• The [( ± 10), 2 ] t laminate exhibited the highest o e stress concentration value and the 
[( ± 45) 12 ], laminate the lowest value. At gauss planes near an interface (Z = H/12 
and Z = 5H/12), the endochronic theory predicted a lower peak <x 9 stress concen- 
tration in the [( ± 10) 12 ], and [( ± 45) 12 ], laminates and a higher peak a 8 stress 
concentration in the [( ± 20) 12 ], and [( ± 30) 12 ], laminates. Both analyses predicted 
the peak a e stress concentration to occur at the hole edge for the [( ± 10) 12 ]„ lami- 
nate and at points a small distance away from the hole edge for the other laminates. 

• The [( ± 30) 12 ] s laminate exhibited the highest x ei stress concentration and the 
[( ± 10) 12 ], laminate the lowest. In all laminates, linear theory predicted significantly 
higher peak t 0z stress concentrations than the endochronic theory. The peak t 0z 
stress concentration occured at the hole edge. 

• The cr z concentration profiles oscillate between positive and negative values, al- 
though the magnitudes are very small. The endochronic theory profiles are very 
similar to the linear theory profiles except near the hole edge where the endochronic 
theory predicts much larger o z stress concentrations than linear theory. 

5. The predicted stress profiles in material principle coordinates were examined at the hole 

edge, for — 90° ^ 0 ^ 90°. The highlights of these profiles are listed below: 

• The linear theory and endochronic theory a, stress profiles are very similar for all 
laminates except the [( ± 45) 12 ], laminate where the endochronic theory predicts a 
40% higher peak a, stress concentration. The endochronic theory predicts higher 
peak ct, stress concentrations than the linear theory for the [( ± 10) 12 ]„ and 
[( ± 45) 12 ], laminates and lower than linear theory concentrations for the 
[( ± 20) 12 ], and [( ± 30) 12 ], laminates. 
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• The linear theory and endochronic theory a 2 stress profiles are similar in the 
[< ± 10)„], and [( dt 20) 12 ], laminates and relatively dissimilar in the [( ± 30) 12 ], and 
[( ± 45) 12 ], laminates. The endochronic theory a 2 stress profiles are much more 
jagged than the linear theory profiles. 

• The linear theory and endochronic theory a 3 (a t ) stress profiles differ significantly for 
ail the laminates except the [( ± 10) 12 ], laminate. The endochronic theory <j 3 stress 
profiles exhibit radical jumps and show stress magnitudes much higher than the lin- 
ear theory profiles. 

• The linear theory and endochronic theory interlaminar shear stress a 4 (t 23 ) profiles 
are similar for all the laminates except the [( ± 20) 12 ] s laminate where the profiles 
show a significant difference in magnitude. 

• The linear theory and endochronic theory interlaminar shear stress <r 5 (x 13 ) profiles 
agree reasonably well for all laminates except the [( ± 20) 12 ], laminate. The 
endochronic theory peak c 5 stress concentrations are significantly lower than the 
linear theory concentrations in all laminates except the [( ± 20) 12 ], laminate. 

• The linear theory and endochronic theory shear stress o, (t 12 ) profiles are similar for 
all the laminates. The endochronic theory <y s stress profiles are much more jagged 
than the smooth linear theory profiles. In all of the laminates, the peak ct s stress 
concentration is lowest for endochronic theory analysis. 

6. Cylindrical coordinate o e and t ei stress profiles were examined at the hole edge for 

-90° £ 9 £ 90°. The highlights of the o e and t 8i stress profiles are listed below. 

• The <j 0 stress profiles are very similar to the a, stress profiles and show the same 
trends. 
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• The t 9i stress profiles show the maximum T e , stress to occur at Q=0° for all laminates. 
The linear theory and endochronic theory stress profiles are similar in shape, but 
differ in peak magnitude. In all laminates, the endochronic theory predicted lower 
peak t 0i stress concentrations than the linear theory. 


6.4 Failure 


A failure analysis of the unnotched and notched angle-ply laminates was performed using the 
tensor polynomial failure criterion. The tensor polynomial was used to predict the fracture 
load and, for the notched laminates, the location and mode of failure. Failure away from the 
hole edge was examined by applying the tensor polynomial at radii greater than the hole ra- 
dius for -90° ^ 0 £ 90°. Finally, the effect of the ratio w/d on failure was investigated for the 
[{ ± 10), 2 ], and [( ± 20) 1Z ], laminates. The important feature of the failure analysis are out- 
lined below. 

1. Theoretical fracture load predictions agreed very well with experiment for ail of the 
unnotched laminates except the [ ± 20], laminate. For the [ ± 20], laminate, theory 
gave non-conservative failure predictions. An in-plane shear strength value of 14 ksi gave 
better fracture load predictions than a shear strength value of 9.4 ksi. The endochronic 
theory predicted o“'* stresses nearly the same as linear theory, but predicted ej 1 * strains 
more accurately than linear theory. 

2. Fracture load predictions for the notched laminates were very conservative. A shear 
strength value of 14 ksi gave better fracture load predictions than a shear strength value 
of 9.4 ksi; however, the fracture load predictions were still very conservative. The 
endochronic theory and linear theory a“'* stresses were very similar. However, the 
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endochronic theory e“'* strains did not correlate with experiment significantly better than 
the linear theory ey" strains because the conservative fracture loads did not allow the 
laminates to be deformed significantly into their softening regions. 

3. The material nonlinear endochronic theory fracture predictions for the notched laminates 
did not correlate with experiment significantly better than the linear theory predictions. 
Therefore, the cause of the conservative fracture loads is likely due to finite element 
stresses near the hole edge that are too high and the application of the tensor polynomial 
at the hole edge. Laminates often exhibit local failure (micro-cracking) near notches be- 
fore fracture; therefore, the tensor polynomial may actually be predicting local failure and 
not fracture. 

4. Application of the tensor polynomial at the hole edge gave very accurate predictions of 
the locations and modes of failure. When a shear strength of 9.4 ksi was used, the 
endochronic theory gave failure modes and locations identical to those of the linear the- 
ory for all laminates. However, when a shear strength of 14 ksi was used, the 
endochronic theory and linear theory predicted different failure locations and different 
tensor polynomial terms for the [( ± 30) 12 ], and [( ± 45) 12 ] s laminates. In the 
[( ± 30 ) 12 ] s laminate, the endochronic theory predicted a different failure mode than the 
linear theory using a shear strength of 14 ksi. 

5. Application of the tensor polynomial at various radii (or distances from the hole edge) for 
-90° ^ 0 £ 90° shows the lowest fracture loads to occur at the hole edge. The profiles 
of far field ay" stress versus the distance D/A do not converge at any D/A value; therefore, 
a characteristic distance does not seem to exist. As D/A increases, the predicted failure 
locations change in both 0 and Z dimensions. The endochronic theory and linear theory 
predict identical failure locations as a function of D/A for all laminates except the 
[( ± 45) 12 ], laminate. 
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6. [( ± 10) 12 ] t and [( ± 20) 12 ] t laminates with w/d ratios of 10 and 5 were tested 1 and ana- 

lyzed. For a decrease in w/d from 10 to 5, experiments revealed a reduction in strength 
of approximately 18.5% and a reduction in ultimate failure strain of approximately 20% 
for both laminates. Both linear theory and endochronic theory gave reductions in strength 
and ultimate failure strain approximately equal to the 11.1% reduction in cross-sectional 
area. 
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Appendix A. Data Reduction 


Appendix A. Data Reduction 


For a given test, a straight line representing a "best fit" through the data is required. The data 
reduction in Chapter 2 uses a least squares fit to best approximate this straight line. A 
straight line is represented by the equation y = mx + b where m is the slope and b is the 
y-intercept on a cartesian coordinate system. If x is designated to be the independent vari- 
able, then for a given x data point, there exists both a y data point (called y-observed or y obs ) 
and a y value corresponding to the equation y = mx + b The difference between these values 
is represented by the equation d = y obl — mx + b . 

It is desired to minimize this difference over the entire range of data points, but since d may 
be a positive or negative quantity, the difference d cannot simply be summed over each point. 
The least squares method sums the square of the difference (function f below) over each point. 
This may be represented in functional form as 

f(m,b) = £d? (A.1) 

i = 1 


where 


d| = y ( - mX| - b 


In order to minimize this function of the two variables m and b, the derivative of the function 
with respect to each variable is set equal to zero and the resulting two equations are solved 
simultaneously. 


JL = JL = o 

dm db 


(A.2) 


df = 
dm 


£ 2(y, - mx, - b)x, = 0 


i = i 



- 2(y ( - mXj - b) = 0 
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Equations (A. 2) may be expanded to yield: 


n n 

m £ X| + bn = £ y ( 
i=i 1=1 


(A. 3) 


n „ n n 

mlXj + b I X| = I x ( y, 
l=i i=i i=i 


In matrix form, equations (A. 3) exhibit the form: 


n 

I X, 

i = 1 

n 9 

I xf 

i = 1 


n 

n 

I Xj 
i = 1 


n 

Svi 

n=i 

, i x .y.| 

i=i 


(A. 4) 


Equations (A.4) yield a "best fit" value for the slope m and the y-intercept b. 
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The following Tables are referenced in Chapter 2 for determining the endochronic constants. 
They contain slopes and y-axis intercepts from various lamina tests. 


Table B.l. Tables for Determining An, v 12 , A 22 , and A 44 

(a) An Determined from 0° Uniaxial Compression Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

An 

1 

19,509,667.9 

13.17 

5.12566E-08 

2 

19,500,513.2 

-48.80 

5.12807E-08 

Average 

19,505,090.5 

-17.82 

5.12687E-08 


(b) v 12 Determined from 0° Uniaxial Compression Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

V 12 

1 

-0.3500225 

-6.81225E-06 

0.3500225 

2 

-0.3159658 

-1.69590E-05 

0.3159658 

Average 

-0.3329942 

-1.18856E-05 

0.3329942 


(c) A 22 Determined from 90° Uniaxial Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

a 22 

1 

1,525,254.8 

150.04 

6.55628E-07 

2 

1,527,561.5 

-55.12 

6.54638E-07 

Average 

1,526,408.2 

47.46 

6.55133E-07 


(d) A 44 Determined from 2-3 Plane losipescu Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

A 44 

1 

472,208.9 

-31.94 

2.11771E-06 

2 

445,017.1 

51.94 

2.24710E-06 

Average 

458,613.0 

10.0 

2.18049E-06 
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Table B.2. Tables for Determining A*si Bn, n,, n 2 , and B° 22 


(a) A«« Determined from 1-2 Plane losipescu Tests and Off-Axis Tests 


losipescu Tests 

Slope m (psi) 

Y-intercept (psi) 

As® 

1 

812,268.9 

1.46 

1.23112E-06 

2 

903,552.3 

0.04 

1.10674E-06 

Off-Axis 

Slope m (psi) 

Y-intercept (psi) 

A 6 6 

10 

839,663.2 

0.00 

1.19095E-06 

15 

887,303.5 

0.00 

1.12701E-06 

30 

789,924.2 

0.00 

1.26594E-06 

45 

836,014.9 

-0.05 

1.19615E-06 

Average 

844,787.2 

0.48 

1.18373E-06 


(b) Bn and ni Determined from 0° Uniaxial Compression Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

n i 

Bn 

1 

2.403743 

-35.79 

2.403743 

2.85762E-16 

2 

2.492074 

-36.89 

2.492074 

9.55283E-17 

Average 

2.447909 

-36.34 

2.447909 

1.90645E-16 


(c) n 2 Determined from 90° Uniaxial Compression Tests 


Test# 

Slope m (psi) 

Y-intercept (psi) 

n 2 

1 

2.84790 

-35.62 

1.8479 

2 

4.07145 

-48.00 

3.0714 

Average 

3.459675 

-41.81 

2.4597 


(d) B° 22 Calculated Using n 2 from 90° Uniaxial Tests 


Test# 

Y-intercept 

B° n x SJf 984 

1 

-41.9532 

2.08449E-18 

2 

-41.7110 

2.65574E-18 

Average 

-41.8321 

2.37010E-18 
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Table B.3. Tables for Determining n 6 and S 66 


(a) n 8A Determined from 1-2 Plane losipescu Tests and Off-Axis Tests Using the Average 
Value of A m 


losipescu Tests 

Slope k 6 

In C 6 

n 6A " 1 

1 

3.84135 

-41.29 

2.84135 

2 

4.04939 

-43.18 

3.04939 

Off-Axis Tests 

Slope k« 

In C« 

"6A = k 6 - 1 

10 

2.99530 

-33.08 

1.99530 

15 

2.96310 

-32.33 

1.96310 

30 

3.73866 

-39.83 

2.73866 

45 

3.68246 

-38.59 

2.68246 

Average 

3.54504 

-31.62 

2.54504 


(b) n 6B Determined from 1-2 Plane losipescu Tests and Off-Axis Tests Using Individual 
Values of A 6 « 


losipescu Tests 

Slope k e 

In C« 

n 6B = lU - 1 

1 

3.84140 

-41.14 

2.84140 

2 

4.18540 

-44.69 

3.18540 

Off-Axis Tests 

Slope ke 

In C 6 

n 6B = ke - 1 

10 

2.99770 

-33.08 

1.99770 

15 

O 07400 
t - . \j i -rdvj 

on co 

A 07400 
1 .CM tOU 

30 

3.61290 

-38.47 

2.61290 

45 

3.58900 

-37.77 

2.58900 

Average 

3.53355 

-37.96 

2.53355 


(c) S«« Calculated Using n 6A for 1-2 Plane losipescu Tests 


Test# 

Y-intercept 

$66 

1 

-38.60 

1.81042E-13 

2 

-38.40 

2.12236E-13 

Average 

-38.50 

1.96639E-13 
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Appendix C. Criterion For Positive Definiteness 


Consider the following matrix: 


S 11 

Si2 

Si2 

0 

0 

0 

s 12 

S 22 

S 23 

0 

0 

0 

S i2 

®23 

s 22 

0 

0 

0 

0 

0 

0 

S44 

0 

0 

0 

0 

0 

0 

S 55 

0 

0 

0 

0 

0 

0 

S ss 


To be positive definite, the above matrix must meet the following criteria: 

8„ > 0 (B.2) 


Sli S 12 

S« S 22 


> 0 


(B.3) 


11 

S-12 

S 12 



•12 

s 22 

®23 

> 0 

(B.4) 

; 12 

CO 

CM 

CO 

S 22 




Expanding equations (B.2) through (B.4) gives the following conditions: 


1. 


8 ,, 


> 0 


2 . 


S„S; 


22 


S? 2 > 0 


3- S^ss* sy + 2sys 23 S 22 ) > 0 
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Appendix D. Finite Element Implementation of 
Inversion Symmetry 


Appendix D. Finite Element Implementation of Inversion Symmetry 


A finite element model is a set of N equations with N unknowns. The inversion symmetry 
conditions introduce M new equations, reducing the number of equations and unknowns to N 
“ M. Applying the inversion symmetry condition to an anisotropic plate in the modeling stage 
reduces the geometry of the model from full to half. Any line passing through the center of 
symmetry which splits the plate into two half sections may be used. The symmetry line cho- 
sen for modeling the plate with a central hole of this study is the y-axis {Figure D.1). 

A simpie example is used here to show how a full finite element model is reduced to a half 
model. Consider an anisotropic plate under compression. A simple four element grid with 4 
node elements represents the finite element model as shown in Figure D.1. For the sake of 
simplicity, assume that each node has only one degree of freedom, that in the y-direction. 


The global set of equations is written 



(D.1) 


where Vj represent the unknown nodal displacements. 

Applying the inversion symmetry relationship to the displacements, using the Y-axis as the 
line of symmetry, 
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Substituting equations (D.2) into (D.1) gives the following set of equations. 



(D.3) 


Equation set (D.3) has only five unknown displacements -- four of the unknown displacements 
were eliminated using inversion symmetry. Rearranging equation set (D.3) yields: 
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Reducing equation set (D.4) to the five independent equations and unknowns yields: 


<*22 + *8®) 

<*23 + *78) 

<*M - K m ) 

<*26 + K 4 „) - (K 24 + Kg,) 

1 

* 

K> 

+ 

1 

i 


\ 

V 2 

<*23 + *78) 

(*33 + *77) 

<*35 - *57) 

(*47 + Kag) 

0 


V 3 

<*25 + *58) 

<*35 ~ *57) 

*55 

<*56 - *45> 

* 

1 

i 

< 

v 5 > 

<*48 + *26) - <*24 + K 69 ) 

<*47 + *36) 

<*56 - *45) 

<*44 + K«) 

<*,4 + *68) 


V. 

""3> 

CO 

* 

1 

+ 

1 

0 

<*59 “ *15) 

<*14 + *68) 

<*11 + *89) 


v V J 



Because the stiffness terms are inversion symmetric as well, equation set (D.5) may be sim- 
plified even further. Note that 


*88 = *22 

*48 = *26 
*12 = *89 


*78 = *23 
*24 = *68 
*77 = *33 
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^57 “ ^35 K47 - 


(D.6) 


^58 = K 2 5 K 56 = K45 

K44 = K 66 K 15 = K59 

•^14 = K 6 9 K„ = Kgg 


Therefore, 


2K 22 

CO 

CM 

* 

CM 

0 

2(K 2 g Kgg) 

r — " — ~ — 

O) 

CO 

CM 

1 

(S 

1 


CO 

CM 

* 

CM 

CO 

CO 

CM 

0 

2K 36 

0 

1 V,j 

1 

CO 

u_ 

CM 

0 

0 

K 55 

0 

0 

V 

- 

0 

2(K 26 — Kgg) 

2K 36 

0 

§ 

CM 

2K 69 

r 6 ! 

1 1 

! °( 

0> 

CO 

CM 

1 

0 

0 

o> 

(O 

* 

CM 

2Kgg 

( v ». 

) 1 

W 


(D.7) 


Equation set (D.7) include stiffness and displacement terms from elements 1 and 2 only. Of 
course, a 2 can be factored from equations 1,2,4, and 5 in equation set (D.7). 

This simple example shows how a full plate finite element model may be reduced to that of 
a half plate model (Figure D.2). In the above example, the line of symmetry was taken as the 
Y-axis. Alternatively, the X-axis, a line through node points 3, 5, and 7, or a line through points 
• 1, 5, and 9 could have been used as the line of symmetry. 

Of course, in modeling the full plate with a half plate mesh, the finite element program must 
be modified internally and the input data modified appropriately. This is discussed next. 
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D.1 Finite Element Development 


Independent and dependent nodes can be defined for the half plate mesh shown in Figure D.2. 
Those nodes associated with the independent degrees of freedom are called independent 
nodes. All other nodes are dependent nodes. In the plate in Figure D.2, node 8 {or node 2) 
is the only dependent node because its displacement is dependent upon the node 2 dis- 
placement (or the node 2 displacement is dependent upon the node 8 displacement). 

In order to implement the inversion symmetry condition, the finite element program must 
know which nodes are independent and which nodes are dependent. This is done in the input 
data. Each node is labeled with an integer valued flag variable -- either a 0, 1, or -1. The 
variable 0 indicates that this node is independent. The variable 1 indicates that the node is 
independent, but has a dependent node which corresponds to it. The variable -1 indicates the 
node is dependent and corresponds to one of the nodes marked with a 1. The dependent 
nodes are matched to their appropriate independent nodes by their global coordinate re- 
lationship. 

The program assigns each node to a global equation number. Dependent nodes are assigned 
to the same global equation number as their corresponding independent nodes. In the as- 
sembly procedure, the stiffness terms of the dependent nodes are added to the stiffness terms 
of the independent nodes. However, before the assembly of the global stiffness matrix, the 
elemental stiffness terms are modified in those elements which contain dependent nodes. 

For elements with dependent nodes, the stiffness and load matrices are modified by transf- 
ormations of the form 
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[k] = [T] T [K][r] (D.8) 

{P} = [T] T {P} 


where [K] and {P> are the element stiffness and load matrices and [rj is the transformation 
matrix. 


For element 2 in the mesh of Figure D.2, the elemental equations would appear as 



(D.9) 


or 


[k 2 ] {V > = {P 2 } 


The transformation matrix would look like 


in 


10 0 0 
0 1 0 0 

0 0-10 
0 0 0 1 


{D. 10] 


After the transformation, the transformed stiffness and load matrices would look like 
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[R] 


K ?1 ^12 

i/2 i^2 

*12 *22 

- k^ 3 — k| 3 
Ku K ^ 4 


— K 2 
*13 

Ku 

CO 
CM CM 

1 

K24 

K33 

-kL 

- 4 . 

K44 


(D.1 1) 



For element 1, the elemental equations would appear as: 



[K 1 ] {V > = {P 1 } 


Assembling the element stiffness and force matrices to form the global stiffness and force 
matrices yields 
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(k!i 

+ ^33) 

*12 

<*« 

- 

*13> 

CO 
CM CM 

1 

r- r- 

— K 2 
*34 

*12 

*22 

*23 

K24 

0 

(*13 

-*13> 

*23 

(*33 

+ 

*1l) 

{K34 + K? a ) 

*14 

<*14 

-K23) 

*24 

(*34 

+ 

*12) 

'cm 

CM CM 
+ 

k | 4 

- 

*34 

0 

*14 

*L 

K44 




Transforming the elemental stiffness terms into global stiffness terms yields 


k 22 

*23 

0 

(*26 — K 68 ) 

*89 

*23 

*33 

*35 

*36 

0 

0 

*35 

*55 

*56 

*59 

(*26 — Kgg) 

*36 

*56 

*66 

*69 

-k 89 

0 

*59 

*69 

*99 



(D. 


These global equations are identical to those of equation set (D.7) in the previous section. 
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To compare scalar operations to vector operations, consider the following subroutine. 


SUBROUTINE XT(N) 

COMMON A(100),B(100),C(100),D{100) 
COMMON X(100),Y(100) 

DO 10 1 = 1,2 

X{l) = A<I)*B(I) + C(I)*D(I-1) 

Y(l) = ABS(X(I)) 

10 CONTINUE 

RETURN 

END 


The sequence of scalar operations in this loop are as follows: 


D(1)*C(2) Register, 

A(2)*B(2) Register, 

Register, + Register, Register K 

Register,* X(2) 

ABS(Register K ) Y(2) 

D(2)*C(3) Register, 

A(3)*B(3) Register 

Register, + Register, Register,* 

Register* X(3) 

ABS(Register K ) Y(3) 


D(N-1)*C(N) Register, 

A(N)*B(N) Register 

Register, + Register Register,* 

Register,* X(N) 

ABS(Register K ) Y(N) 


The Cyber 205 vector operation sequence is shown below: 


D(1)*C(2) TEMP,(1) 

D(2)*C(3) TEMP,(2) 
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D(N-1)’C(N) 


TEMP,(N-1) 


A(2)*B(2) TEMP 2 (1) 

A(3)*B(3) TEMPj(2) 

* 


A(N)*B(N) TEMP 2 (N-1) 

TEMP^I) + TEMP 2 (1) X(2) 

TEMP 1 (2) + TEMP 2 (2) X(3) 

* 


TEMP,(N-1) + TEMP 2 (N-1) — X(N) 

ABS(X(2)) Y(2) 

ABS(X{3)) Y(3) 

* 


ABS(X(N)) Y(N) 


As the example above shows, the vector processor performs its operations on vectors (or 
arrays) such as D(N-1)*C(N). It also may perform the above operations at the same time, re- 
ferred to as parallel processing. Obviously, this processing is much speedier than the scalar 
processing. 


However, not all do-loops can be vectorized. An array which is a function of previous values 
of itself may not be vectorized, as shown in the example below. 


SUBROUTINE XT(N) 

COMMON A(100),B(100),C(100),D(100) 
COMMON X(100),Y{100) 

DO 10 1 = 2, N 

X(l) = A(I)*B(I) + C(I)*X(I-1) 

Y(l) = ABS(X(I)) 

10 CONTINUE 

RETURN 

END 
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The above operation may only be computed one step at a time because the value of X(l) de- 
pends upon X(l-1). Vector processing would give erroneous results in this case. Cyber 205 
vector processing restrictions are: (1) Arguments must be determinable prior to the operation, 
(2) Operands must be addressable under certain conventions, and (3) Array elements must 
be contiguous. 
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